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1  Introduction 


The  use  of  adhesives  in  bonded  joint  concepts  has  given  rise  to  a  myriad  of  analytical  and  numerical  treat¬ 
ments  to  ascertain  the  state  of  stress  in  critical  regions  of  the  bond  to  predict  ultimate  joint  strength  and 
service  life.  Analytical  approaches  have  developed  since  the  classic  work  of  Golan  and  Reissner  [1]  in  which 
the  elastic  response  of  a  single-lap  joint  is  analyzed  using  simplifying  assumptions  for  the  stress  state  and 
strain-displacement  relations.  Subsequent  efforts  have  involved  developing  analyses  accounting  for  various 
joint  parameters  such  as  nonisotropic  adherends  and  bondline  thickness  effects  [2-4],  adhesive  flexibility 
[3,5],  geometric  and  material  nonlinearity  of  the  adhesive  and  adherends  [6-10],  effect  of  adhesive  spew  fillet 
[11,12],  stress  singularities  [13-18],  environmental  effects  [19],  combined  loading  [20,21],  and  various  joint 
configurations  such  as  double  lap,  tapered  and  cylindrical  [16,22,23]. 

For  all  the  success  of  analytical  approaches  in  elucidating  the  effect  of  joint  parameters  on  joint  behavior, 
an  inescapable  limitation  is  the  complexity  of  incorporating  all  significant  parameters  in  a  single  analytical 
framework  for  the  analysis  of  general  joint  configurations.  Such  restrictions  have  limited  the  applicability  of 
specialized  analytic  approaches  for  the  design  of  practical  bonded  joint  concepts.  A  preferred  approach  is 
a  finite  element  based  methodology  which  can  provide  a  robust  numerical  technique  to  model  any  bonded 
joint  configuration  while  theoretically  capable  of  accounting  for  all  joint  variables. 

Various  finite  element  formulations  specialized  for  the  analysis  of  adhesive  joints  have  been  reported  in  the 
literature.  A  simple  shear-spring  element  was  developed  in  [24]  which  neglected  the  inclusion  of  normal 
stresses  and  was  used  to  determine  stress  intensity  factors  in  cracked  and  uncracked  adherends.  A  spe¬ 
cialized  element  was  developed  in  [25]  yielding  a  combined  normal  and  shear-spring  representation  of  the 
adhesive.  Carpenter  [26]  formulated  1-D  and  2-D  elements  to  be  joined  to  plate  element  approximations 
of  the  adherends -which  incorporate  parameters  allowing  various  simplifying  assumptions  regarding  adhesive 
behavior  to  be  represented  in  the  element  formulations  in  order  to  numerically  simulate  analytical  results. 
An  extension  of  Carpenter’s  approach  to  3-D  applications  is  presented  in  [27]  but  the  element  derived  is  a 
rod-type  formulation  not  suited  for  3-D  elastic  continuum  representation.  A  mixed  finite  element  is  discussed 
in  [28]  which  explicitly  enforces  interface  stress  conditions  but  details  of  the  formulation  are  not  presented 
and  solutions  are  limited  to  butt  joint  configurations  without  comparison  to  reference  solutions.  In  the 
present  analysis,  the  aim  has  been  to  fully  develop  a  finite  element  basis  for  a  versatile  numerical  approach 
wherein  the  complete  equations  of  2-D  and  3-D  elasticity  are  used  to  accurately  represent  the  elastic  con¬ 
tinuum  and  obtain  results  for  actual  adhesive  joint  configurations  with  realistic  boundary  conditions  and 
applied  loading.  Of  critical  importance  in  the  analysis  of  adhesive  joints  from  a  design  standpoint  are  the 
accurate  prediction  of  bondline  stresses  at  the  ends  of  the  joint  which  tend  to  be  maximum  and,  hence,  the 
initiation  sites  for  adhesive  failure.  The  numerical  prediction  of  stresses  in  these  regions  is  complicated  by 
the  presence  of  boundary  layer  singularities  at  the  joint  ends  due  to  the  free-edge  and  material  mismatch 
effects  and  by  the  nonlinear  geometric  and/or  nonlinear  material  behavior  of  the  joint  under  normal  design 
load  conditions.  The  present  study  has  been  conducted  to  evaluate  the  application  of  the  hybrid  stress 
element  method  [29-33]  in  special  layered  element  configurations  to  accurately  obtain  bondline  stresses.  The 
selection  of  hybrid  stress  element  formulations  is  motivated  by  the  versatility  inherent  in  the  method  by 
allowing  independent  assumptions  regarding  the  stress  and  displacement  field  to  be  incorporated  into  ele¬ 
ment  designs.  The  independence  of  stresses  may  be  used  in  layered  element  configurations  to  strictly  enforce 
stress  continuity  at  the  layer  interfaces  which  is  difficult  in  displacement-based  element  formulations.  An 
additional  capability  is  the  modelling  of  zero  traction  boundary  conditions  which  may  be  used  to  accurately 
represent  true  stress  states  at  the  bond  ends.  The  present  efTort  is  thus  aimed  at  fully  exploring  various 
modifications  to  the  hybrid  stress  variational  basis  while  formulating  elements  with  various  node  configu¬ 
rations,  layer  combinations,  order  of  interpolants  and  order  of  stress  expansions  to  determine  the  optimum 
element  design  for  application  to  predicting  adhesive  joint  stresses.  Numerous  multi-layered  2-D  and  3-D 
elements  formulations  are  detailed  herein  to  assess  their  relative  performance  in  the  prediction  of  interface 
stresses  in  adhesive  joint  configurations.  With  the  emphasis  on  optimum  element  formulation,  the  present 
analysis  is  limited  to  linear  elastic  material  behavior  and  a  single  common  adhesive  joint  configuration  is 
selected  to  benchmark  element  performance.  Towards  this  aim,  a  single-lap  joint  configuration  is  selected 
for  detailed  comparative  analysis.  Additional  enhancements  such  as  the  representation  of  nonlinear  material 
and  geometric  behavior  in  the  presence  of  singular  stress  fields  are  suggested  for  future  development.  The 
performance  of  the  developed  adhesive  finite  elements  are  compared  with  reference  solutions  obtained  using  a 
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highly  discretized  model  of  common  single-lap  adhesive  joint  configurations  using  higher-order  displacement- 
based  elements  and  demonstrate  greater  efficiency  and  accuracy  of  the  developed  hybrid  stress  elements.  All 
specialized  ‘adhesive  elements’  were  run  using  the  ABAQUS  commercial  finite  element  program  through  the 
development  of  a  user-defined  subroutine.  The  complete  source  code  and  a  description  of  the  use  of  these 
elements  is  presented  in  [34]. 

2  Modelling  Stresses  in  Adhesive  Joints 

In  the  joining  of  structural  components,  mechanical  fastening  techniques  have  the  undesirable  effect  of  cre¬ 
ating  highly  complex  stress  fields  and  excessive  stress  intensification  in  the  region  of  the  fastener.  As  an 
alternative,  adhesive  bonding  has  been  increasingly  applied  as  a  joining  technique  to  alleviate  the  stress 
concentration  by  distributing  the  load  transfer  region.  Although  eliminating  the  concentration  of  binding 
forces  in  localized  areas,  however,  adhesive  bonding  gives  rise  to  a  highly  complex  material  system  in  the 
adhesive  layer.  In  the  local  region  of  the  adhesive-adherend  interface  exists  a  critical  interphasal  region  in 
which  various  chemical  species  formed  in  the  reaction  of  adhesive  and  adherend  materials,  together  with 
possible  reagents  used  as  surface  treatments  to  enhance  bonding,  create  a  complex  stratification  of  material 
properties.  Characterizing  the  mechanical  properties  of  these  transitional  zones  is  important  to  quantify 
possible  reductions  in  stiffness  ’and  strength  compared  to  the  bulk  properties  of  the  adhesive  and  adherend 
material.  The  modelling  of  adhesive  layers  thereby  presents  the  common  problem  of  scale  in  representing 
physical  structures:  at  a  microscopic  level  the  representation  of  local  stresses  in  the  complex  material  system 
of  the  interface  is  of  interest  from  the  standpoint  of  fracture  and  failure;  at  a  higher  level,  the  global  response 
of  an  entire  adhesively  bonded  component  is  needed  to  determine  overall  stress  distribution  under  applied 
service  loads.  It  is  the  prediction  of  macroscopic  behavior  which  is  of  interest  in  the  current  study;  global 
stress  predictions  are  required  to  provide  an  estimate  of  local  stresses  to  be  used  in  microstructural  models 
for  thedetailed  analysis  of  complex  micromechanical  response  or  in  ultimate  strength  criteria  for  predicting 
the  initiation  of  joint  failure.  The  analysis  of  a  structure  or  component  containing  an  adhesive  layer  thus 
requires  that  the  variation  in  the  microstructure  of  the  joint  be  neglected  and  gross  mechanical  properties 
used  for  a  numerical  representation.  The  material  property  variation  in  the  local  interface  region  are  replaced 
by  the  assumption  of  a  homogeneous  continuum  in  order  to  predict  the  transfer  of  normal  and  shear  loads 
between  the  adherends  through  the  adhesive  layer.  The  critical  stresses  are  thus  generated  as  normal  and 
tangential  components  along  the  adhesive/adherend  interface,  normally  in  the  region  of  the  ends  of  the  joint. 
At  a  structural  level  of  analysis,  complications  arise  in  the  finite  element  analysis  of  adhesive  stresses  due,  in 
general,  to  the  inherent  mathematical  idealizations  made  in  representing  an  adhesive  joint  and,  in  particular, 
to  the  assumption  of  linear  elastic  material  response.  The  interaction  of  free  edges  and  material  mismatch 
at  the  bond  ends  tend  to  give  rise  to  stress  singularities  and  may  be  resolved  through  the  incorporation  of 
assumed  singular  stress  fields  in  the  element  formulations  or,  most  simply,  the  region  in  the  area  of  the  sin¬ 
gularity  may  be  regarded  as  nonconvergent  and  material  yielding  is  assumed  a  priori.  Nonlinear  viscoelastic 
and  plastic  material  behavior  of  the  adhesive  under  service  loads  commonly  tend  to  reduce  peak  stresses  at 
the  joint  ends  through  redistribution  of  joint  loads.  A  linear  analysis  of  adhesive  stresses  using  appropriate 
failure  prediction  methods  such  as  the  average  stress  criteria  is,  therefore,  conservative  in  predicting  joint 
strength. 

A  numerical  simulation  of  critical  joint  stresses  requires  an  accurate  modelling  of  adhesively  bonded  sections. 
This,  in  turn,  requires  an  accurate  representation  of  the  elastic  continuum  of  the  adhesive  and  adherends 
along  the  bond  line.  Standard  finite  element  formulations  cannot  simultaneously  model  the  discontinuous 
jump  in  material  properties  while  maintaining  continuity  of  stresses  across  the  adherend/adhesive  interface. 
Nonlinear  geometric  and/or  material  behavior  may  be  accounted  for  through  a  variety  of  solution  proce¬ 
dures  accounting  for  large  displacements  and  viscoelastic/plastic  material  models.  The  focus  of  the  current 
research  is  to  assess  the  application  of  the  hybrid  stress  technique  through  the  development  of  novel  element 
configurations  to  accurately  account  for  material  property  discontinuity  while  maintaining  stress  continu¬ 
ity  at  the  critical  adhesive/adherend  interfaces  in  2-D  and  3-D  bonded  joint  configurations.  The  primary 
goal  is  to  rigorously  assess  possible  special  element  designs  involving  different  orders  of  assumed  displace¬ 
ments,  node  configuration,  number  of  layers,  assumed  orders  of  stress  expansions  and  constraints  imposed 
on  the  stress  field  to  determine  the  optimum  element  formulation.  In  order  to  limit  the  amount  of  joint 
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variables  to  discern  comparative  element  behavior,  analyses  are  restricted  to  a  single  common  adhesive  joint 
configuration  incorporating  linear  elastic  material  behavior  and  small  displacement  response  to  applied  loads. 


3  Variational  Basis  for  Specialized  Adhesive  Elements 

The  hybrid  stress  finite  element  method  is  used  in  the  present  study  to  develop  specialized  layered  ele¬ 
ment  formulations.  The  structure  of  the  constituent  element  matrices  are  defined  by  the  Hellinger-Reissner 
functional  in  which  stresses  and  displacements  are  assumed  as  independent  field  quantities.  Alternative 
statements  of  the  basic  hybrid  stress  variational  basis  may  be  utilized  to  optimize  element  performance 
through  assumptions  made  in  the  selection  and  treatment  of  independence  of  assumed  stress  fields.  Two 
basic  versions  of  the  variational  basis  are  investigated  and  are  referred  to  as  the  ‘full’  and  ‘partial’  hybrid 
stress  method  which  are  detailed  below. 


3.1  Full  Hybrid  Stress  Method 

The  complete  Hellinger-Reissner  functional  may  be  stated  as 

Hr  =  /[(  —  l/2)<rrS<r  +  <rT(Lu)]dv  —  f  uTtds  (1) 

Jv  J$ 

where  <r  is  the  assumed  stress  field,  S  is  the  material  compliance  matrix,  u  is  the  assumed  displacement 
field,  L  is  the  differential  operator  relating  strains  to  displacements,  and  t  are  applied  surface  tractions  over 
a  portion  of  the  element  boundary,  s. 

The  assumed  stresses  for  the  i,h  layer  may  be  represented  by 


<r'  =  P'0‘  (2) 

where  P‘  is  a  matrix  of  polynomial  terms  and  /3‘  is  a  vector  of  undetermined  expansion  coefficients. 

The  displacement  field  is  assumed  over  the  ith  layer  domain  as 

u’  =  N‘q'  (3) 

where  N*  are  compatible  isoparametric  displacement  shape  functions  and  q'  are  nodal  displacements.  Ne¬ 
glecting  applied  tractions  and  substituting  (2)  and  (3)  into  (1)  yields 


or 


where 


Hfl  =  /  [(-l/2)/J<TP<TSiP</Ji  +0<TP<T(LN)‘]q<]dt' 

J  v 

n’R  =  (— l/2)/3,rH</3'  +  /3'T  G'q’ 

H‘  =  /  P‘T  S’P'dti 
Jv 

G*  =  f  P'r(LN’)dv  =  f  P'T  B’di; 

Jv  Jv 

Performing  local  assembly  over  the  N  element  layers  yields 

N 

n a  =  £  IT*  =  (-  1/2)/3tH/3  +  ^TGq 

i  =  l 

Seeking  a  stationary  value  of  the  functional  by  taking  the  first  variation  with  respect  to  /3  yields 

0  =  H-IGq 

By  substituting  0  into  (8),  the  variation  with  respect  to  q  yields  the  element  stiffness  matrix  as 

K  =  GtH~*G 


(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

(10) 
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3.2  Partial  Hybrid  Stress  Method 

An  augmentation  of  the  full  hybrid  stress  variational  basis  can  be  made  in  which  not  all  stress  components 
are  represented  by  independently  assumed  stress  expansions.  The  resulting  basis,  called  the  partial  hybrid 
stress  method,  incorporates  independently  assumed  stresses  for  selected  components  while  the  remaining 
components  are  obtained  directly  from  the  assumed  displacement  functions  as  in  classical  displacement- 
based  element  formulations.  This  gives  the  total  stress  field  as  the  sum 


<T  -  <TX  +  <Ti  (11) 

Assumed  independent  stress  components  for  the  ith  layer  are  given  by 

<ri=P*/3‘  (12) 


while  the  remaining  stress  components  are  derived  from  assumed  displacements  in  the  iitl  layer  and  are 
expressed  as 

<r‘2  =  C’B'q*  .  (13) 

where  C*  is  the  material  stiffness  matrix  partition  relating  dependent  stresses  to  strains,  B‘  is  the  strain- 
displacement  matrix  and  q‘  is  the  vector  of  nodal  displacements.  The  variational  statement  of  the  partial 
Hellinger-Reissner  functional  neglecting  applied  surface  tractions  may  be  stated  as 


n*R  =  / [(l/2)qiTB<TC*B*q‘+^,TPiTB*q*  -  ( l/2)/3'T PiT S'P'p^dv 
Jv 


(14) 


where  the  constituent  matrices  for  G*  and  H‘  are  defined  above.  Performing  local  assembly  over  the  element 
layers  yields 


N  r 

n„  =  J2  n*fl  =  /  [( l/2)qTBTCBq  +  /3TGq  -  (l/2)/3TH d\dv 
»=1  Ju 


(15) 


Taking  variations  with  respect  to  the  independent  field  quantities,  the  contributions  to  the  element  stiffness 
matrix  due  to  assumed  stresses  is  given  by 


Kl=GrH-‘G  (16) 

and  the  contribution  due  to  assumed  displacements  is  given  by 

K2  = 

The  total  element  stiffness  matrix  is  therefore  given  by  the  sum 

K  =  K!+K2  (18) 

This  formulation  minimizes  the  computational  cost  associated  with  the  hybrid  method  and  offers  the  benefit 
of  increasing  the  accuracy  of  selected  stress  components  while  allowing  other  components  to  be  computed 
less  accurately  from  assumed  displacements.  The  partial  hybrid  formulation  is  used  to  assess  2-D  element 
performance  for  potential  application  in  3-D  layered  element  formulations  in  which  the  large  number  of  in¬ 
dependent  stress  modes  used  for  maintaining  complete  expansions  may  be  minimized  while  suppressing  zero 
energy  modes.  An  important  aspect  of  the  partial  variational  statement  is  that  it  precludes  the  enforcement 
of  all  field  equilibrium  conditions  on  the  element  stresses  which  may  limit  element  performance  in  specific 
applications. 

Additional  variations  to  both  of  the  above  hybrid  variational  basis  may  be  made  through  the  selective 
relaxation  of  layer  equilibrium  and  interface  stress  continuity  constraints. 
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BTCBdv 


(17) 
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4  Layered  Element  Formulations 


In  using  the  hybrid  approach,  the  constituent  matrices  of  individual  layers  are  locally  assembled  and  stress 
continuity  conditions  may  be  explicitly  enforced  at  layer  interfaces  resulting  in  novel  element  formulations 
ideally  capable  of  representing  the  stress  state  at  material  interfaces.  In  the  historical  development-  of  the 
hybrid  stress  method,  the  improvement  in  element  behavior  has  been  offset  by  an  ostensibly  greater  com¬ 
putational  cost  in  forming  element  stiffness  coefficients.  This  has  resulted  in  the  convention  of  utilizing 
minimum  expansions  for  assumed  stress  components  to  minimize  additional  ccet  over  displacement-based 
element  formulations.  Although  the  basic  assumption  that  the  hybrid  stress  technique  inherently  requires 
a  greater  computational  cost  over  displacement-based  element  formulations  has  been  disputed  recently  in 
References  [35,36],  the  development  of  specialized  elements  presupposes  a  limited  use  of  such  elements  in 
a  complete  structural  finite  element  model.  Thus,  the  computational  efficiency  of  specialized  elements  are 
regarded  as  of  secondary  importance  and  complete  emphasis  is  placed  on  developing  element  formulations 
which  provide  the  most  accurate  stress  prediction.  An  additional  aspect  of  the  hybrid  method  has  been  the 
problem  of  selecting  assumed  stress  fields  which  suppress  spurious  zero  energy  modes  in  the  resulting  element 
stiffness  matrix.  This  has  been  solved  by  the  development  of  rational  procedures  for  selecting  stress  expan¬ 
sions  based  on  stresses  assumed  in  mapped  or  natural  coordinates  [32,33].  In  this  approach  a  contravariant 
transformation  based  on  Jacobians  computed  at  the  element  centroid  is  used  to  map  natural  stresses  into 
physical  coordinates.  For  the  elements  developed  in  the  current  study,  a  minor  geometric  restriction  is 
imposed  which  requires  that  all  element  layers  are  rectangular.  This  requirement  simplifies  the  mapping 
function  and  allows  stresses  to  be  assumed  in  physical  coordinates.  In  the  selection  of  assumed  stresses, 
complete  stress  expansions  for  each  layer  of  various  orders  are  used  to  assess  element  characteristics  and 
performance. 

The  various  treatments  of  selected  stress  fields  and  the  performance  of  2-D  and  3-D  specialized  layered 
elements  -  or  ‘adhesive  elements’  -  are  presented  beiow. 

4.1  2-D  Layered  Element  Configurations 

Figure  1  depicts  typical  layered  element  configurations  used  in  the  present  study  to  represent  a  2-D  elastic 
continuum.  Nodal  variables  are  restricted  to  u  and  v  translational  degrees  of  freedom.  The  isoparametric 
displacement  fields  presented  below  are  assumed  independently  within  each  element  layer;  superscripts  de¬ 
noting  layer  designation  have  been  dropped  for  clarity. 
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Figure  1.  Typical  2-D  layered  element  configurations. 
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As  shown  in  Figure  2,  a  local  coordinate  system  is  assumed  within  each  layer  with  the  origin  situated 
at  the  layer  centroid.  For  the  various  element  configurations  shown  in  Figure  1,  the  layer  displacement  shape 
functions  are  presented  below. 


Figure  2.  Local  and  mapped  layer  coordinate  systems  in  2-D  and  3-D  elements. 

For  a  linear  element  with  4  nodes  per  layer,  the  isoparametric  displacement  functions  are  given  by 


(19) 


The  mapping  between  physical  and  natural  coordinates  is  given  by 


*  =  <*0  +  +  Os’?  +  a3^f] 

y  =  I>q  +  +  b^t]  +  63^^ 


where 


«o  bo 

1  1  11' 

*1  y  1 

«i  bi 

_  1 

-1  1-11 

*2  y-i 

0-2 

~  4 

-1-1  11 

*3  ito 

O3  63 

1-1-11 

.  x*  y*  . 

(20) 


(21) 


For  a  higher-order  element  with  8  nodes  per  layer,  the  isoparametric  displacement  functions  u9  are  given  by 


(22) 


where 


N»  =  \ 


’  Nl  1 

(i  -  0(1  -*/)(-! -€-*») 

N2 

2(1  —  (*)(1  —  v) 

(1+0(1 -«?)(-! +^-1) 

/V4 

2(1  —  0(i  -  n2) 

iVs 

’  —  ' 

2(1+  0(1  -r,2) 

N6 

(1 -0(1 +  i?)(-l -*  +  >?) 

n7 

2(1  -  OKI  +  0 

Ns 

(1  +  0(1 +  »»)(-! +€  +  »?) 

(23) 


The  mapping  between  physical  and  natural  coordinates  is  given  by 

x  =  a0  +  al(  +  a2rf  +  a3(r/+a4(2  +  asr/2  +  a6£3r/  +  aT(r/2 
y  =  60  +  6i£  +  bit)  +  63^7  +  b^2  +  b$t)2  +  bet2Ji  + 


(24) 
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where 
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For  a  mixed-order  element  with  displacements  assumed  as  a  linear  function  of  the  r  coordinate  and  a  quadratic 
function  of  y  resulting  in  6  nodes  per  layer,  the  isoparametric  displacement  functions  u,  are  given  by 


(26) 


where 


Ni=  < 


The  mapping  between  physical  and  natural  coordinates  is  given  by 


JVi 

iKi-O(i-i) 

Ni 

1(1  +0(i-  i) 

2(1  —  0(1  —  I2) 
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►  —  < 
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‘  - - - 

x  =  a0  +  ai4  +  <»2>?  +  «3£»7+a4»?2  +  «5£»?2 
y  =  b0  +  bit  +  bit)  +  63£q  +  64jj2  +  6s£l2 
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Stress  equilibrium  conditions  for  the  ith  layer  are  given  by 


6x  6y  6x 

Continuity  conditions  for  stresses  at  the  layer  interfaces  are  given  by 


=0 


(27) 


(28) 
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xv'n=h. 


- 


=  <r\ 


yy  4 


_  r<+i| 


(30) 


(31) 


*»  'n=-h, 


In  the  use  of  the  partial  hybrid  formulation,  only  the  normal  and  shear  stress  components  are  approximated 
through  independent  expansions.  A  priori  enforcement  of  layer  equilibrium  is  reduced  to  the  condition 


Sr‘  6ff' 

_£]£.  + —Slit  =  0 

ox  by 


(32) 


while  the  interface  continuity  conditions  remain  those  given  in  (31). 
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4.2  3-D  Layered  Element  Configurations 

Representative  3-D  element  configurations  incorporating  multiple  layers  are  depicted  in  Figure  3. 


Figure  3.  Typical  3-D  layered  element  configurations. 


The  isoparametric  displacement  functions  u,  for  each  layer  are  given  by 


u,  ]  ,  a  Ui 

|  =|E(1+^)(1  +  ^)(1+C*C)|  w*  |=N.q  (33) 


The  isoparametric  mapping  between  physical  and  natural  coordinates  is  given  by 


where 


x  =  a0  +  arf+a2rj  +  a3(  +  arfij  +  a5£C  +  +  0"£tC 

y  —  b0  +  brf  +  63^  +  63C  +  brfrj  -f  brf£  +  fr6*?C  +  b?£i]C 

Z  =  C0  +  erf  +  C2T)  +  C3(  +  crff)  +  crf£  +  C6T)C  +  crfljC 
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07  67  C7  -1  1—1  1  1  —  1  1  —  1  _  *8  2/8  *8 


(34) 


(35) 
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Stress  equilibrium  conditions  for  the  i,h  layer  are  given  by 


60'-,  St^  St*. 

— ££.  +  _£!  +  _££.  =  0 
Sx  Sy  6z 

St*  So*  Sr* 

■  **  j - VL  j - LL  =  A 

Sx  Sy  6z 

St*.  St*  Bo*.. 

1 f+-#+-sf  =  0 


(36) 


Continuity  conditions  of  stresses  at  the  layer  interfaces  are  given  by 


-  t<=_fc1+, 

Ti*i<  =  k, 

=  T6*M(=_*,+| 

(37) 

In  the  partial  hybrid  approach,  only  the  ot,  rix  and  ryi  components  are  based  on  independently  assumed 
stress  expansions.  The  enforcement  of  equilibrium  throughout  the  i,k  layer  domain  is  limited  to 


6t'x,  8rh 

Sx  6y  6z 


=  0 


(38) 


The  continuity  conditions  of  stresses  at  the  layer  interfaces  remain  the  same  as  given  in  (37). 


5  Benchmark  Solutions  for  Adhesive  Stresses 


A  variety  of  2-D  and  3-D  layered  element  formulations  have  been  derived  to  assess  the  effect  of  element 
layers  and  order  of  stress  expansion  on  the  accuracy  of  stress  prediction.  All  elements  have  been  supported 
for  use  in  the  ABAQUS  commercial  finite  element  program  through  a  user-defined  subroutine  [34].  The 
developed  elements  are  benchmarked  using  a  2-D  single-lap  joint  shown  in  Figure  4  and  two  3-D  single-lap 
joints  with  a  rectangular  and  tapered  pianform  configuration  as  shown  in  Figures  5  and  6.  In  order  to  assess 
the  accuracy  of  stress  prediction  using  specialized  elements,  a  reference  solution  was  obtained  by  using  a 
highly  discretized  2-D  finite  element  model  incorporating  a  large  number  of  higher-order  displacement- based 
8-node  plane-strain  elements.  Regular  meshes  were  used  for  each  model  domain  with  the  discretization 
indicated  using  the  notation  (nr, ny)  to  indicate  the  number  of  elements  used  in  the  x  and  jr  directions, 
respectively.  The  boundary  conditions  and  applied  loading  used  in  both  2-D  and  3-D  single-lap  models 
are  depicted  in  Figure  7.  The  roller  supports  at  the  adherend  ends  are  applied  over  the  first  few  elements 
to  enforce  a  zero  curvature  condition.  The  uniform  applied  load,  P,  is  distributed  as  a  tensile  load  over 
the  joint  end  in  the  upper  adherend  while  fixity  conditions  are  applied  at  the  opposite  end  of  the  lower 
adherend.  Although  numerous  analytical  approaches  have  been  presented  in  the  literature  for  the  analysis 
of  single-lap  joints,  all  incorporate  simplifying  assumptions  in  order  to  yield  closed-form  solutions  which 
are  incompatible  with  the  approximation  framework  of  finite  element  theory.  As  an  example,  governing 
differential  equations  describing  a  single-lap  joint  are  presented  in  Reference  [15]  and  were  solved  to  generate 
an  analytical  solution  for  bondline  stresses.  Assumptions  made  in  deriving  the  governing  equations  include 
zero  transverse  deformation  of  the  adherends,  constant  stresses  through  the  adhesive  thickness  and  joint 
flexibility  primarily  due  to  deformation  in  the  adhesive  layer.  As  shown  in  Figures  8  and  9,  the  analytical 
solution  does  not  match  the  solution  obtained  numerically  using  a  finite  element  based  approach  -  most 
notably  in  the  region  near  the  singular  point  present  in  the  finite  element  analysis  assuming  linear  elastic 
material  response  -  and  hence  is  unsuitable  for  the  present  study. 
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Figure  5.  Geometry  of  a  3-D  single-lap  joint  model  with  rectangular  planform. 


Figure  6.  Geometry  of  a  3-D  single-lap  joint  model  with  a  tapered  planform. 
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Figure  7.  Boundary  conditions  used  to  define  the  single-lap  joint  model  supports. 


ao  1.0  2.0  3.0  4.0  5.0  6.0 
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Figure  8.  Analytic  and  finite  element  prediction  of  c ryy  stresses  along  bondline  in  a  single-lap  joint. 


DISTANCE  ALONG  BONDUNE.  X 

Figure  9.  Analytic  and  finite  element  prediction  of  rxy  stresses  along  bondline  in  a  single-lap  joint. 


Id  addition  to  providing  the  reference  solution  for  developed  2-D  elements,  the  plane-strain  solution  is  used 
directly  as  a  reference  solution  for  the  3-D  single-lap  joint  with  rectangular  planform  assuming  negligible 
variation  of  stresses  through  the  joint  width.  For  the  tapered  3-D  single-lap  joint,  the  bondline  stresses  in 
the  reference  solution  are  scaled  according  to 


<r„  =  (1  +  x/l) 

*z,  =  (1  +  *//)**» 

where  x  is  the  axial  coordinate  and  1  is  to  total  bond  length. 


(39) 


Numerous  element  formulations  and  their  performance  are  detailed  in  the  following  sections.  Different 
combinations  of  element  layers,  node  configuration,  element  order,  stress  field  selection  and  the  use  of  ‘full’ 
and  ‘partial’  hybrid  stress  variational  bases  are  assessed.  Due  to  the  significantly  greater  degree  of  develop¬ 
mental  effort  in  3-D  element  formulations,  basic  element  performance  is  assessed  in  2-D  configurations.  The 
optimum  element  design  is  then  carried  over  to  the  development  of  a  3-D  element.  The  development  of  a 
generic  3-D  element  with  optimum  stress  prediction  capablilities  is  the  ultimate  aim  of  the  current  research 
to  provide  a  versatile  numerical  tool  for  the  modelling  of  any  arbitrary  3-D  joint  configuration. 


6  2-D  Special  Adhesive  Elements 


Various  2-D  elements  are  formulated  and  their  accuracy  in  predicting  joint  stresses  is  assessed  through  a 
relative  comparison  with  the  reference  solution.  As  explained  above,  a  detailed  study  of  2-D  element  configu¬ 
rations  is  performed  tc  ascertain  the  optimum  element  configuration  based  on  node  configurations,  number  of 
element  layers,  order  of  assumed  displacement  shape  functions,  order  of  assumed  stress  expansions  and  vari¬ 
ations  in  stress  field  constraints  for  the  development  of  an  optimum  3-D  element  formulation  which  involves 
a  significantly  greater  degree  of  complexity  in  basic  formulation.  The  various  2-D  element  formulations  and 
their  accuracy  in  predicting  joint  stresses  in  a  single-lap  configuration  are  discussed  in  detail  below.  In  the 
discussions  of  particular  element  performance,  statements  regarding  the  hybrid  stress  finite  element  method 
are  presented  which  may  provide  general  insight  into  the  use  of  this  variational  basis  for  specialized  layered 
finite  element  formulations  for  other  intended  applications.  To  assess  the  convergence  of  the  derived  ele¬ 
ments,  four  different  models  utilizing  10,  25,  50  and  100  elements  along  the  bondline,  were  constructed  and 
compared  to  the  reference  solution.  These  models  were  specifically  constructed  with  nonoptimal  meshes  in 
the  use  of  available  degrees  of  freedom  in  the  local  region  of  the  adhesive  bond  to  assess  element  performance 
in  models  with  high  aspect  ratios  and  significant  jumps  in  the  degree  of  mesh  refinement  between  regions  in 
the  model.  The  models  are  depicted  in  Figure  10  with  the  discretization  indicated  in  various  regions  by  the 
notation  (NtxNr)  where  Ni  denotes  the  number  of  elements  used  along  the  coordinate  axis.  The  width 
of  the  joint  is  assumed  to  be  unity  and  isotropic  material  properties  are  selected  as: 

Adherend  :  E  =  69000.0  n  —  0.32 

Adhesive  :  E  =  3000.0  fi  =  0.36 

All  stresses  are  normalized  as  er'j  =  Oij/ffrej  where  <rrtj  =  P/A  in  which  P  is  a  uniformly  applied  tensile 
load  and  A  is  thfe  cross-sectional  area  of  the  adherend  end. 
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N.r  2  FOR  2-LAYERED  HYBRID  ELEMENTS 
i  1  FOR  3-LAYERED  HYBRID  ELEMENTS 


Figure  10.  Discretization  used  in  the  2-D  models  to  assess  solution  convergence. 


6.1  The  H2L6N  Element 

The  H2L6N  element  is  formulated  as  a  6-node,  two  layer  element  as  shown  in  Figure  11.  This  element 
is  constructed  specifically  to  provide  a  more  accurate  representation  of  stresses  at  material  interfaces  by 
explicitly  enforcing  stress  continuity  of  the  <ryy  and  rry  stress  components.  In  modelling  the  adhesive  layer, 
two  H2L6N  elements  are  used  through  the  adhesive  thickness  and  are  joined  to  displacement-based  4-node 
plane  strain  elements  used  to  model  the  adherends.  The  use  of  the  H2L6N  element  in  modelling  an  adhesive 
layer  is  depicted  in  Figure  12.  4-node  plane-strain  displacement-based  elements  -  designated  CPE4  -  were 
obtained  from  the  ABAQUS  element  library  and  used  to  model  the  adherends. 
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Figure  11.  The  H2L6N  element. 


2-LAYER  ADHESIVE 
ELEMENTS 


Figure  12.  The  use  of  special  two-layer  elements  in  modelling  adhesive  bonds. 

Three  versions  of  H2L6N  were  developed  to  assess  the  order  of  stress  expansions  on  element  performance.  In 
all  formulations  field  equilibrium  is  enforced  a  priori  together  with  the  interface  continuity  conditions.  The 
number  of  independent  strain  modes,  nf,  may  be  computed  as 

nc  ”  ™do/  ftrbm  (40) 
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where  nj0j  are  the  number  of  element  degrees  of  freedom  and  nr*m  are  the  number  of  rigid  body  modes.  For 
the  H2L6N  element,  the  bilinear  isoparametric  displacement  held  yields  nt  =  12  -  3  or  9  independent  strain 
modes.  Version  I  incorporates  complete  linear  expansions  for  each  layer  resulting  in  a  stress  field  possessing 
10  independent  stress  modes  which  is  sufficient  to  suppress  any  zero  energy  modes  in  the  resultant  element 
stiffness  matrix.  The  stress  expansions  are  given  by: 


°lx 

= 

0\  +  02*  +  03V 

_1 

= 

04  +  0s*  -  07y 

riy 

= 

06  +  07X  -  07y 

<*lx 

= 

03  +  09*  +  0\OV 

= 

04  +  0sx  —  07(h\  +  h2  +  y) 

Try 

= 

06  —  0102  +  07*  —  09(02  +  y) 

Version  II  incorporates  complete  quadratic  expansions  for  the  assumed  stresses  yielding  18  independent  stress 
modes  given  by 

alx  =  Pi-  fc*  +  0ay  -  204 *y  -  0s*2 /2  +  06y 2 

ffjy  =  07  +  20a[(hi  +  h2)x  -  xy]  +  fax  +  fcy{hi  -y/2)  -  0ioy  +  0nh?y  +  0i2*2 

=  0\3  +  0s*{y~  0\)  +  0iOx  -  f3nh2z  +  fay  +  fox2  +  fay2  (42) 

<rl*  =  0\4  —  0is*  +  06y  —  20tfiy  -  flux2 /2  +  fay2 

<*yy  —  07  +  0sh\/2  —  0io(.h\  +  h2  +  y)  +  0n[hili2  +  (ho  —  y‘)/2]  +  09*  —  2/?8xy  +  0\ox2 

Txy  =  013  + 02^1  +  04^1  +  0is(l*2  +  y)  +  0l~(y2  —  h?)  +  010*  +  0ll*y  +  08*2 

Version  III  incorporates  complete  cubic  expansions  for  the  assumed  stresses  with  28  independent  modes 
given  by 

<f\s  =  0i-  *02  +  V03  -  2xy04  -  x20s/2  +  y2 06  -  Zxy207  -  yx20a  -  x3/?9/3  +  y30io 
ffyy  =  (— hf/2  +  hxy  -  y~/2)0$  +  (h\y  —  2h3/3  -  y3/3)/J8  +  (hi  +  h2  -  y)0i3  +  0\2  + 
h2(y  —  hj  —  hn/2)0i4  +  h?(hi  +  h2/ 3  —  y)0\s  +  x(2hxy  —  h\  —  y~)0 9  +  x0\6  + 

2x(h\  +  h2  —  y)0n  +  xh2(2y  —  2/»i  —  ho)0  is  +  3x"(hi  +  h2  —  y)0 19  +  x~0xx  +  x302o 
Txy  =  (y  —  hl)02  +  (y2  ~  h2)04  +  (y3  ~  h,3)0T  +  021  —  hn022  +  02023  —  00024  +  x(y  —  h\)0s  +  (43) 

*(y2  ~  02)0s  +  X013  -  h2x0XA  +  h\x0 15  +  x2(y  -  hi  )09  +  x20l7  -  hox20xa  +  x30x9 
°lt  =  02S  -  X022  +  y02 6  ~  2 xy023  ~  X2/?h/2  +  y2 027  ~  3 xy2024  -  yx20i 5  -  x30ia/i  +  y302S 
o\y  =  0\2  +  *016  -  y0 13  -  2 xy0n  +  x20xx  -  y2/ 20lA  -  xy20 xa  -  3yx2/?19  +  x3/?20  -  y30 )5/3 
Txy  =  021  +  *013  + y022  +  *y0\4  +  *~017  +  y2023  +  *y20lS  +  y*20is  +  *3019  +  y3024 

The  element  behavior  in  the  prediction  of  <ryy  and  rzy  stresses  utilizing  the  various  complete  orders  of  as¬ 
sumed  stress  expansions  are  depicted  below.  The  four  basic  models  using  10,  25,  50  and  100  II2L6N  elements 
along  the  bondline  are  compared  with  the  reference  solution  to  assess  convergence  characteristics. 

Figures  13  and  14  show  the  distribution  of  <ry  and  rry  stresses  across  the  bondline  using  the  II2L6N  el¬ 
ement  incorporating  linear  stress  expansions.  Rapid  convergence  to  the  reference  solution  is  observed  for  all 
but  the  highly  coarse  10-element  model. 
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Figure  13.  H2L6N  Version  I  prediction  of  <ryy  stresses  along  the  bondline. 


Figure  14.  H2L6N  Version  I  prediction  of  rxy  stresses  along  the  bondline. 


Figures  15  and  16  show  the  distribution  of  tryy  and  rty  stresses  across  the  bondline  using  the  H2L6N 
element  incorporating  complete  quadratic  stress  expansions. 
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Figure  15.  H2L6N  Version  II  predictions  of  ayy  stresses  along  the  bondline. 


Figure  16.  H2L6N  Version  II  prediction  of  rxy  stresses  along  the  bondline. 


Rapid  convergence  to  the  reference  solution  is  again  observed  with  improved  convergence  behavior  demon¬ 
strated  in  the  coarse  mesh  10-element  model. 

Figures  17  and  18  show  the  distribution  of  <ryy  and  rxy  stresses  across  the  bondline  using  the  H2L6N 
element  incorporating  complete  cubic  stress  expansions. 
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Figure  18.  H2L6N  Version  III  prediction  of  rfy  stresses  along  the  bondline. 


Comparisons  of  the  various  models  show  that  all  element  formulations  yield  excellent  predictions  for  both  the 
normal  and  shear  stress  distributions.  The  cubic  order  expansions  used  in  version  III  provide  little  improve¬ 
ment  over  the  quadratic  fields  incorporated  in  version  II  which  indicates  the  diminishing  return  in  adopting 
stress  fields  of  increasingly  higher  order.  Also  demonstrated  is  the  improvement  in  the  coarse  10-element 
model  with  increasing  order  of  assumed  expansions  for  stresses. 
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In  order  to  quantify  the  improvement  in  the  coarse-mesh  element  behavior  in  the  10-element  model,  a 
comparison  of  the  eigenvalues  of  the  hybrid  element  stiffness  matrices  using  complete  linear,  quadratic  and 
cubic  stress  expansions  is  performed  using  the  element  geometry  shown  in  Figure  19.  A  comparison  is  also 
made  to  a  displacement-based  formulation  in  which  two  separate  4-node  elements  where  locally  assembled 
to  compare  with  the  layered  hybrid  elements.  The  element  geometry  is  selected  to  show  the  effect  of  a  high 
aspect  ratio  layer  which  may  be  encountered  in  representing  the  adherend/  adhesive  layers  in  a  coarse  model. 
The  material  properties  used  were  E  —  1000  and  \i  =  0.25  and  the  joint  width  was  taken  as  unity.  The 
results  depicted  in  Table  1  show  that  by  adopting  a  higher-order  stress  representation,  ‘weak’  deformation 
modes  present  in  the  linear  field  element  are  removed  which  can  be  related  to  the  improvement  in  coarse- 
mesh  behavior. 


Figure  19.  Sample  configuration  of  H2L6N  containing  a  high  aspect  ratio  layer. 


Table  1.  Comparison  of  eigenvalues  obtained  in  the  H2L6N  element  incorporating  complete 
linear,  quadratic  and  cubic  stress  expansions  with  a  displacement-based  element. 


A 

Quadratic  Expansions 

1  »JH»!f!THi7SITHSS35!HI 

1 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0' 

0.0 

0.0 

0.0 

0.0 

0.0 

2.389470 

59.08495 

167.7284 

324.0885 

6.608349 

218.5657 

257.2027 

346.4270 

250.7947 

428.2339 

469.4921 

644.1404 

355.4981 

464.6144 

545.8709 

651.8773 

8 

536.3158 

713.2400 

853.5973 

1219.328 

9 

642.8778 

892.4405 

1039.350 

13791.84 

10 

905.2318 

930.4139 

1125.610 

40050.97 

11 

1233.321 

1961.596 

2418.629 

40289.39 

12 

1924.831 

2324.951 

2675.814 

120303.3 

The  conclusion  from  the  above  study  is  that  the  minimum  stress  expansion  used  in  version  I  is  adequate  to 
suppress  zero  energy  modes  yet  is  deficient  as  compared  to  higher  orders  of  stress  expansions  in  terms  of 
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coarse  mesh  accuracy.  The  reason  for  the  significant  change  in  element  characteristics  with  increasing  order 
of  assumed  stress  fields  is  due  to  the  constraints  imposed  in  formulating  a  layered  element.  These  conditions 
constrain  the  element  stress  field  such  that  the  inner  product  of  assumed  functions  comprising  the  stress 
space  with  those  defining  the  strain  space  is  highly  dependent  on  the  order  of  stress  expansions  used.  The 
increase  in  the  spanning  of  the  strain  space  is  measured  by  the  eigencharacteristics  of  the  resulting  element 
stiffness  matrices  and  represents  the  elastic  strain  energy  associated  with  particular  element  deformation 
modes.  Another  observation  is  that  the  extremely  ‘stiff  deformation  modes  evidenced  in  the  displacement- 
based  ‘layered’  element  formulation  due  to  the  high  aspect  ratio  of  the  bottom  element  layer  is  not  present  in 
the  hybrid  element  formulations.  The  comparison  of  element  performance  based  on  complete  quadratic  and 
cubic  stress  field  expansions  -  versions  II  and  III  -  show  negligible  difference  indicating  that  the  quadratic 
stress  field  is  optimum  for  this  element  configuration. 

To  compare  the  H2L6N  element  predictions  with  conventional  finite  elements,  figures  20  and  21  show  results 
using  standard  displacement-based  elements.  The  elements  used  are  the  CPE4  plane-strain  quadrilaterals 
available  in  the  ABAQUS  library.  In  the  displacement-based  model,  two  4-node  plane  strain  elements  were 
used  to  model  the  adhesive  layer  through  the  thickness  utilizing  the  same  number  of  degrees  of  freedom  as  in 
the  layered  hybrid  elements.  Because  the  displacement-based  solution  is  incapable  of  representing  continuity 
of  stress  components  across  the  interface,  stresses  were  extrapolated  from  Gauss  points  and  averaged  at 
nodes  for  all  connected  elements.  As  shown,  the  solution  for  bond  interface  stresses  actually  converges  away 
from  the  reference  solution  with  increasing  mesh  refinement  along  the  bond  line.  The  improvement  in  stress 
prediction  using  the  layered  hybrid  stress  approach  is  related  to  the  ability  of  explicitly  imposing  the  stress 
continuity  conditions  at  the  interface  and  equilibrium  constraints  with  each  layer  domain. 


DISTANCE  ALONG  BONDLINE,  X 

Figure  20.  CPE4  prediction  of  <ryy  stress  distribution  along  the  bondline. 
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Figure  21.  CPE4  prediction  of  rxy  stress  distribution  along  the  bondline. 

Stress  recovery  at  points  away  from  the  critical  bond  interface  region  show  a  closer  agreement  between  the 
hybrid  and  displacement-based  elements.  Figures  22  through  25  show  stresses  recovered  at  element  centroids 
and  indicates  that  the  special  layered  adhesive  elements  and  displacement-based  elements  converge  to  the 
same  solution  within  material  domains  while  the  hybrid  elements  demonstrate  greater  efficiency  and  accu¬ 
racy  in  stress  prediction  along  critical  domain  boundaries  where  material  properties  are  discontinuous  such 
as  along  the  adhesive/adherend  interface. 


Figure  22.  H2L6N  prediction  of  ayy  stress  distribution  computed  at  element  centroid. 
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Figure  23.  H2L6N  prediction  of  rzy  stress  distribution  computed  at  element  centroid. 


Figure  24.  CPE4  prediction  of  <ryy  stress  distribution  computed  at  element  centroid. 
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Figure  25.  CPE4  prediction  of  Txy  stress  distribution  computed  at  element  centroid. 


In  the  hybrid  stress  method,  a  limitation  principle  exists  in  which,  for  elements  with  rectangular  geome¬ 
try,  the  element  stiffness  coefficients  become  identical  to  those  obtained  from  a  purely  displacement-based 
approach  if  no  constraints  are  imposed  on  the  assumed  stress  fields.  An  examination  was,  therefore,  un¬ 
dertaken  to  qualify  the  improvement  in  hybrid  element  performance  over  displacement-based  elements  by 
assessing  the  the  relative  effect  of  selectively  enforcing  layer  equilibrium  and  interface  continuity  conditions 
in  element  formulations.  Figures  26  and  27  show  the  prediction  of  <ryy  and  Tly  stresses  along  the  bondline 
using  complete  quadratic  stress  field  expansions  in  which  only  field  equilibrium  constraints  were  enforced  in 
each  layer.  The  resulting  stress  field  contains  24  independent  stress  modes  given  by 


<7*r  =  01  +  X02  +  y/?3  +  xyfi 4  +  *20s/2  +  y2p 6 

<rxyy  =  Ih  +  *0s  +  y09  +  xyfiio  +  *20u  -  y20s/2 
Txy  =  012  -  *09  -y02  +  *2/05  ~  *201O  ~  2002 

(t\x  =  013  +  *014  +  2/015  +  *2/016  +  *2017/2  +  y'018  (44) 

ffyy  =  /?19  +  *020  +  y021  +  *y022  +  *2023  -  y  017/2 

Tly  =  024  _  *021  —  2/014  +  *y0l7  -  *2  022  ~  !/2014 
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Figure  26.  H2L6N  ayy  bondline  stress  prediction  with  only  equilibrium  constraints  enforced. 


Figure  27.  II2L6N  rxy  bondline  stress  prediction  with  only  equilibrium  constraints  enforced. 


As  shown,  the  normal  peel  stresses,  ayy,  are  accurately  predicted  while  the  shear  stress  distribution  shows  a 
convergence  away  from  the  reference  solution.  Although  not  as  pronounced  as  in  the  use  of  CPE4  elements, 
this  departure  from  the  reference  solution  demonstrates  the  importance  of  enforcing  the  continuity  condi¬ 
tions  at  the  layer  interface  and  qualifies  the  improved  accuracy  obtained  in  the  hybrid  stresss  technique  over 
displacement-based  formulations  when  equilibrium  is  enforced  pointwise. 

Figures  28  and  29  show  the  stress  prediction  of  the  1I2L6N  element  in  which  complete  quadratic  stress 
fields  were  constrained  only  to  satisfy  interface  continuity  conditions.  The  resulting  stress  field  contains  30 
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independent  stress  modes  given  by 

*1,  =  06  +  y20i  +  *J/?J  +  xy03  +  y0A  + 

=  A  +  *(y  -  *i)/*7  +  V09  +  *JAo  -  1  +  *£12  +  y20iz 

ri,  =  0i9  +  y10i4  +  zy0i6  +  x70u  +  y017  +  z0lB  (45) 

°r*  =  070  +  I/V21  +  X2  022  +  Xy023  +  y024  +  */?2S 

=  09  +  (y,-h2)026  +  (h2  +  y)027  +  h\08  +  X20io  +  Xy0n  +  X0i2  +  h20i3 
r*»  —  +  (v3  —  hl)07»  +  X(y  +  Aj)/3 29  +  (/«2  +  y)0M  +  ^20n  +  h\x0n  -f  x20it  +  h\0n  +  x0i& 


DISTANCE  ALONG  BONDLINE.  X 

Figure  28.  H2L6N  <ryy  bondline  stress  prediction  with  only  interface  constraints  enforced. 


DISTANCE  ALONG  BON  DUNE,  X 


Figure  29.  H2L6N  rzy  bondline  stress  prediction  with  only  interface  constraints  enforced. 

As  shown,  both  the  normal  stress  and  the  shear  stress  distribution  is  shown  to  converge  away  from  the 
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reference  solution.  The  conclusion  from  the  above  results  indicate  that  the  imposition  of  both  layer  equilib¬ 
rium  and  interface  stress  continuity  constraints  is  necessary  for  obtaining  optimal  element  performance  in 
predicting  bondiine  stresses. 

The  partial  hybrid  stress  approach  is  next  employed  to  assess  the  effect  of  fully  enforcing  interface  con¬ 
tinuity  of  stress  components  while  incorporating  only  partial  satisfaction  of  the  layer  domain  equilibrium 
conditions.  Two  stress  fields  were  developed  for  <ryy  and  rxy  based  on  different  orders  of  complete  stress 
expansion.  The  normal  stress  component,  a„,  was  calculated  from  the  assumed  displacement  field.  Lin¬ 
ear  stress  expansions  are  sufficient  to  suppress  any  spurious  kinematic  modes  and  the  resulting  stress  field 
contains  6  independent  modes  given  by 

=  07  +  (hi  +  hj  -  y)0i  +  x03 

Tlt  —  0*  +  (if —  ~  (46) 

0f,  =  07  +  -  y0i 

Tzt  —  /%  +  +  y06 

Complete  quadratic  expansions  were  also  developed  containing  12  independent  stress  modes  given  by 

0ft  =  04  +  (hf  —  2hiy+  y3)0i  +  (y  —  hi  —  1*7)02  +  h3(2hi  +  —  2  y)0$  + 

x0s  +  2ar(h,  +  h2  -  y)0o  +  x207 

Tlt  —  0io  +  (y  —  f*i)0a  +  (y2  —  l*i)09  —  f*20n  +  h^0i2  +  (47) 

2x(h\  -  y)0i  -  X02  +  2h7x03  +  x206 
0ft  =  04  +  x0s  +  y02-  2xy06  +  x207  +  y20 3 
r2,  =  0io  -  x02  +  y0i  1  -  2 xy03  +  x206  +  y20 12 

Figures  30  through  33  depict  stress  predictions  using  the  complete  linear  and  quadratic  stress  fields  for  the 
two  independent  stress  components  <ryy  and  rry  as  presented  above. 


DISTANCE  ALONG  BONDLINE,  X 


Figure  30.  H2L6N  <ryy  bondline  distribution  using  partial  hybrid  stress  formulation 
with  linear  stress  fields. 
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Figure  31.  H2L6N  rxy  bondline  distribution  using  partial  hybrid  stress  formulation 
with  linear  stress  fields. 


Figure  32.  II2L6N  <ryy  bondline  distribution  using  partial  hybrid  stress  formulation 
with  quadratic  stress  fields. 


27 


Figure  33.  H2L6N  rx,  bondline  distribution  using  partial  hybrid  stress  formulation 
with  quadratic  stress  fields. 

From  the  above  distributions  it  is  evident  that  the  normal  stress  distribution  is  predicted  with  good  ac¬ 
curacy.  The  results  for  the  shear  stress  distribution  converge  away  from  the  reference  solution  indicating 
the  deficiency  of  the  partial  stress  method  in  accurately  predicting  all  independent  stress  components.  The 
increase  in  the  order  of  expansions  from  linear  to  quadratic  again  demonstrate  the  improvement  in  coarse 
mesh  accuracy  as  shown  by  the  behavior  of  the  10-element  model. 


In  an  actual  bonded  joint,  the  ends  of  the  adhesive  layer  are  traction  free.  Using  conventional  finite  el¬ 
ement  models  gives  rise  to  the  prediction  of  nonzero  shear  and  inpiane  tractions  in  this  region.  In  the  hybrid 
stress  method,  the  ability  to  strictly  enforce  zero  traction  conditions  of  the  rzy  and  oxx  components  at  the 
free  edge  is  obtained  through  a  priori  imposing  the  folbwing  constaints  in  the  assumed  stress  field 


ri„  =0 


=  0 


(48) 


where  a  is  the  element  half-width.  Four  versions  of  an  assumed  quadratic  stress  field  in  which  zero  traction 
is  specified  along  a  particular  face  as  depicted  in  Figure  34  are  developed.  The  full  hybrid  variational  basis 
is  utilized  in  the  element  formulation.  All  stress  fields  possess  15  independent  stress  modes  resulting  in  ‘end’ 
elements  of  correct  rank. 


Figure  34.  H2L6N  element  showing  side  designatbns  on  which  zero  traction  conditions 
are  explicitly  enforced. 
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The  resulting  stress  field  expansions  in  which  zero  traction  conditions  are  enforced  on  various  element  sides 
are  given  below. 

For  zero  traction  on  element  side  FI: 


*1*  =  0i  -  *02  +  v/h-  2*V04  -  *J0s/2  +  y*0t 

ayy  —  08  +  (— h2  +  2hiy  -  y2)/Js/2  +  (y  —  hi  —  h2)07  +  h2(y  —  hi)/?g  + 

2z(hi  —  y)0io  +  x0n  -  h2*/?i2  +  *Vi3  (49) 

ri»  =  (y  -  hi)0 j  +  (y2  -  h\)0A  -  (a  +  x)01  +  (x2  -  a2)0io  - 
h2(a  +  *)09  +  *(y  -  Ai  )0s 
°l,  -  -«(*  +  o)/39  +  (x2  -  a2)0 14  +  y(x  +  a)/315 
avv  =  0s  +  *0n  +  y07  +  xy/?i2  +  *20i3 
TxV  =  -(*  +  a)07  +  (*2-a2)0io  +  y(*  +  a)09 
For  zero  traction  on  element  side  F2:' 

<rlIX  =  — a(x  +  0)09  +  ( x 2  -  a2)0n  +  y(x  +  a)/3i5 
=  08  +  *011  +  y07  +  *y012  +  X2  0\3 

rzy  =  -(*  +  «)07  +  (*2  -  <*2)01O  +  y(x  +  a)09  (50) 

a\s  -  01  -  *02  +  y03  -  2xy/?4  -  z20$/2  +  y206 

ayy  ~  0s  ~  (hj  +  2h2y  +  y")0s/2  + (hi  +  A2  +  y)07  —  hi(h2  +  y)09 — 

2x(h2  +  y)0io  4-  *0n  +  hi*/?i2  +  *20i3 
r«,  =  (h2  +  y)02  +  (y2  -  h^)#,  -  (a  +  *)07  +  (*2  -  a2)^io  + 
hi(a  +  x)0g  +  *(h2  +  y)0s 

For  zero  traction  on  element  side  F3: 

a\x  -  01  ~  *02  +  y03  -  2xy/34  -  x2/?5/2  +  y2/36 

‘’’yy  =  08  —  (h2  —  2A2y  +  y2)05/2  +  (y  —  h2  —  ^2)07  +  h2(y  —  A2)^9 + 

2*(h2  —  y)0io  +  *0n  —  h2*0i2  +  *20i3  (51) 

riy  =  (y  ~  hi)^2  +  (y2  —  h2)/?4  +  (a  —  *)07  +  (*2  —  a2)0io  + 
h2(a  -  *)09  +  *(— hi  +  y)05 
<r2r  =  a(x-a)/?9  +  (x2-a2)/(?M  +  y(x-a)/?i5 
O'yy  =  08 +  *011  +  y07  +  *y012  +  *2013 
r2y  =  (a  -  *)07  +  (*2  -  o2)0io  +  y(*  -  a)/?9 

For  zero  traction  on  element  side  F4: 

°**  =  <>(*  -  o)09  +  (*2  -  a2/0u  +  y(*  -  a)0i5 
ffyy  =  08  +  *01 1  +  y0 7  +  *y012  +  *2013 

riy  =  (a  -  *)07  +  (*2  -  o2)/?io  +  y(*  -  a)A,  (52) 

°\z  =  01  ~  *02  +  y03-  2xy0A  -  x20$/2  +  y20e 

a\y  —  08  —  (h2  +  2h2y  +  y2)/?5/2  +  (hi+h2  +  y)/?7  +  hj(—h2  —  y)09  — 

2*(2h2  +  y)0io  +  *0n  +  h\x0n  +  *20i3 
rzy  =  (h2  +  y)02  +  (y2  -  h\)0 4  +  (a  -  *)/?7  +  (x2  -  a2) 010  + 
hi(*  -  <i)09  +  *(h2  +  y)0s 
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behavior  is,  however,  regarded  as  a  particular  effect  of  the  models  being  used  in  the  present  study;  the  singular 
nature  of  the  normal  stress  solution  is  theoretically  unaffected  through  the  imposition  of  zero  tractions  of 
the  shear  stresses  in  linear  elastic  analysis. 

6.2  The  H2L10N  Element 

The  H2L10N  element  is  configured  as  a  2-layered  element  incorporating  higher-order  displacement  fields  in 
the  bond  thickness  direction.  The  H2L10N  element  configuration  is  depicted  in  Figure  37.  The  interest  in 
this  element  formulation  is  to  assess  the  increase  in  the  order  of  assumed  displacement  fields  selectively  in  the 
direction  of  higher  stress  gradients.  Two  versions  of  this  element  were  derived  utilizing  complete  quadratic 
(version  I)  and  complete  cubic  (version  II)  assumed  stress  fields.  The  assumed  stress  fields  are  identical  to 
the  quadratic  and  cubic  expansions  used  in  II2L6N  and  presented  in  equations  (42)  and  (43).  The  behavior 
of  these  element  formulations  are  presented  in  Figures  38  through  41. 


Figure  37.  H2L10N  element  configuration. 


Figure  38.  H2L10N  Version  I  <ryy  stress  distribution  along  the  bondline. 
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NORMALIZED 


Figure  39.  H2L10N  Version  I  rxy  stress  distribution  along  the  bondline. 


Figure  40.  H2L10N  Version  II  ayy  stress  distribution  along  the  bondline. 
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Figure  41.  I12L10N  Version  II  Txy  stress  distribution  along  the  bondline. 

The  behavior  of  these  elements  demonstrates  two  interesting  effects.  Regardless  of  the  higher-order  strain 
representation  in  the  adhesive  thickness  direction,  the  normal  stress  predictions  are  inferior  to  the  lower- 
order  H2L6N  element  in  predicting  maximum  stresses  in  the  singular  region  of  the  joint.  The  shear  stress 
recovery  in  the  quadratic  stress  version  shows  an  inherent  tendency  of  the  element  to  approximate  a  zero 
traction  condition  at  the  end  of  the  bond.  Increasing  the  order  of  stress  expansion  tends  to  diminish  this 
behavior  and  increases  the  maximum  normal  and  shear  stress  prediction  in  the  region  of  the  singular  point. 
Although,  once  again,  increasing  the  order  of  assumed  stress  expansions  improves  the  results  compared  to 
the  reference  solution,  the  overall  behavior  of  this  element  is  inferior  to  the  linear-order  H2L6N  element. 


6.3  The  H3L8N  Element 

The  H3L8N  element  is  formulated  as  an  8-node,  three  layer  element  as  shown  in  Figure  42.  This  element  is 
formulated  to  model  both  top  and  bottom  interfaces  of  the  adhesive  bond  with  the  inner  layer  representing 
the  adhesive  domain.  The  use  of  the  II3L8N  element  in  adhesive  joint  modelling  is  depicted  in  Figure  43.  In 
anticipation  of  the  high  degree  of  constraint  imposed  on  the  stess  field  to  satisfy  equilibrium  and  continuity 
conditions  in  three  layer  domains  and  two  interfaces  within  a  single  element  configuration,  two  different 
assumed  stress  expansions  of  increasing  order  were  developed.  The  two  versions  of  H3L8N  were  formulated 
incorporating  complete  quadratic  and  cubic  stress  fields  to  assess  the  effect  on  element  performance.  All 
assumed  stress  expansions  provide  sufficient  independent  stress  modes  to  span  the  13  basic  element  strain 
modes  thereby  suppressing  spurious, kinematic  modes.  Stress  continuity  of  the  <ryy  and  rry  stress  components 
is  strictly  enforced  at  each  layer  interface  and  equilibrium  is  enforced  within  each  element  layer. 
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Figure  42.  H3L8N  element  configuration. 


Figure  43.  The  use  of  special  3-layered  elements  in  modeling  adhesive  bonds. 


Version  I  incorporates  complete  quadratic  expansions  for  each  layer  in  which  equilibrium  is  enforced  a  priori 
together  with  the  continuity  conditions.  The  resulting  stress  field  possesses  24  independent  stress  modes 
which  is  sufficient  to  suppress  any  zero  energy  modes  in  the  resultant  element  stiffness  matrix.  The  stress 
expansions  are  given  by 
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a\x  =  0i  +  02x  4-  03V  +  0<xy  -0ix2/2  4-  06y2 

<r\y  =  09  +  0iox  -  0nz(hi  +  h2  -  y)  +  0x2x2  -  03(h\ -2hiy  +  y2)/2  + 

07(y  -hi  —  h2)  +  0»h2(y  -hi  -  h2/2) 

Txy  =  0\3  +  0i(hi  -  y)  4-  04(hi  -  y2)/2  4-  0i4h2  -  0ish$/2  -  0sx(hi  -  y)  - 
07 X  -  08h2x  -  0ux2/2 

o\x  =  016  +  014X  4-  0ny  +  0is xy  -  08x2/2  +  018y2 

avv  ~  09  +  0iqx  +  07y  +  0nxy  +  0x2x2  —  08y2  )2  (53) 

Tlv  =  013  ~  07X  -  014V  +  0nxy  -  0iiz2/2  -  0isy3/2 
<r3x  —  0i9  +  02ox  +  02iy  +  022xy  —  023i2 /2  +  024y2 

°yy  =  09  +  07(h2  +  h3  +  y)  —  0%h2(h3  +  h2/2  +  y)  —  033^3 +  2h2y  +  y2)/2  + 

0iqx  +  0ux(h2  4-  h3  +  y)  +  0X2x2 

T*y  =  013  —  0i4h2  —  0i8h\/2  —  02o(h3  4-  y)  4-  022(hl  -  y2)/2  -  0jx  +  0sh2x  4- 
023x(h3  +  y)  -  0nx2 /2 

Version  II  is  formulated  with  complete  cubic  stresses  to  assess  the  effect  of  higher-order  stress  expansions  on 
element  performance.  The  resulting  stress  field  contains  38  independent  stress  modes  and  is  given  by 

fflx  =  0i+  0i x  +  03y  +  04xy  -  0$x2 / 2  4-  08y2  -  307 xy2  -  foyx2  -  09x3/3  4-  0ioy3 

ffyy  =  0 21  —  03(^i  ~  2hxy  +  y2)/2  —  0g(2h3  -  3 h\y  +  l^)/3  +  023(hx  +  h2  —  y)  —  0\sh2(hi  —  h2/ 2  4-  y)  + 

0iah2(hx  +  h2/3  —  y)  —  09  x(h\  —  2h\y  +  y2)  +  022x  4-  2024x{hx  +  h2  —  y)  — 

0i9x[2h2y  —  h2(2hi  4-  /*2 )]  4-  3028x2(h\  4-  h2  —  y)  4-  02*x2  4-  027X3 
Tly  —  023  +  02 (hi  —  y)  4-  04 (A2  —  y2)/2  —  07(h3  —  y3)  +  0i2h2  —  0\4h\/2  —  0nh\  —  0sx(hx  -  y)  — 

0%z{h\  -  y2)  +  023x  -  0nh2x  +  0iShlx  -  09x2(hi  -y)  +  0nx2  -  0i9h2x2  +  026x3 
=  0n  +  0nx  4-  0i3y  +  014  xy  —  0i$x2 /2  4-  0i6y2  —  30nxy2  —  0i8yx2  —  0X9x3/3  4-  02oVZ 

<*yy  —  02i  +  022X  —  023y  —  2024iy  +  02sx2  —  0i$y2 /2  -  0i9xy2  —  3026yx2  +  027X3  —  0isy3 /3  (54) 

rxy  —  028  +  023X  —  0\2y  +  0i$xy  4-  024X2  -  0i4y2 /2  +  0igiy 2  +  0i9yx 2  +  02qx3  +  0ny3 

fl**  =  029  +  03QX  +  031  y  +  032xy  —  033X2 /2  +  03 J\y2  —  3033  Xy2  —  03$yx2  —  037 X3 / 3  +  033  y3 

ffyy  =  02i  —  023(y  +  h2  +  A3)  —  0ish2(h3  +  Aj/2  —  y)  —  0ish2(h3  4-  Ao/3  +  y)  —  033(h3  +  2hsy  +  y2)/ 2  + 
036(2h3  +  3A|y  —  y3)/3  —  2024x(h2  +  h3  +  y)  —  20i9h2x(h3  +  h2/2  +  y)  + 

037x(h\  —  2A2A3  —  2h2y  —  y2)  +  022x  +  02$x~  —  302ex2(h2  +  A3  +  y)  -t-  027x 3 
Txy  —  ~0i2h2  —  0i4h\/2  4-  0nh3  4-  02a  —  /?3o(A3  4-  y)  4-  032{h\  —  y2)/ 2  4-  /?3s(A|  4-  y3)  4-  023x  4-  0\*>h2x  4- 

0\8h\x  4-  033x(h3  4-  y)  —  036x(h^  —  y2)  4-  024X2  4-  0\9h2x2  4-  037X2{h3  +  y)  +  026X3 


Figures  44  through  47  show  the  stress  predictions  of  interface  bond  stresses  for  both  versions  of  the  H3L8N 
element. 
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jure  44.  H3L8N  Version  I  cryy  stress  distribution  along  the  bondline. 
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stress  field  required  to  obtain  convergent  element  behavior  must  be  contrasted  with  the  use  of  H2L6N  which 
yields  greater  accuracy  while  requiring  only  lower-order  stress  field  expansions. 

An  H3L8N  element  formulation  explicitly  enforcing  zero  traction  conditions  was  developed  and  stress  pre¬ 
dictions  presented  in  Figures  48  and  49.  A  quadratic  stress  field  was  selected  which  contains  21  independent 
modes  as  given  below. 

For  zero  tractions  enforced  on  element  side  FI 

o’L  =  0i  -  +  yfo  -  2xy/?4  -  x20s  +  y206 

=  0s  +  (— hj  +  2yAj  —  y2)0 5/2  +  (y  —  h\  -  h2)07  +  A2(y  —  hx)09  + 

2x(hi  -  y)0io  +  x0n  -  xh20i2  +  x20i3 
Tty  =  (y~  1*1)02  +  (y2  ~  h\)P*  ~(a  +  x)07  +  (x2  -a2)0lo  - 
h2(a  +  x)0g  +  x(y  —  hx)03 
alx  -  ”«(*  +  <*)09  +  (x2  -  a2)0u  +  y(x  +  a)0is 

=  0s  +  x0u  +  y07  +  xy0i2  +  x20i3  (55) 

Txy  =  -(*  +  a)07  +  ( x 2  -  a2)0 10  +  y(x  +  a)09 

fflx  =  0i6-x0i7  +  y0ia-xy0i9-x202O/2  +  y2021 
Oyy  —  0s  —  (h3  +  2Aay  -|-  y2)02o/2  +  (h2  +  h3  +  y)07  —  h2(h3  +  y)09  — 

2  x(h3  +  y)0 10  +  x0n  +  h2x0\2  +  x20x3 

Txy  =  ( h3  +  y)0n  +  ( y 2  -  hl)0i9  -  (a  +  x)07  +  (x2  ~  a2)0iO  +  h2(a  +  x)09  + 
x(h3  +  y)02o 


For  zero  tractions  enforced  on  element  side  F2 

o\x  =  01 -  X02  +V03-  2xy0A  -  x20s  +  y2 0q 

Vyy  =  0a  +  (—hi  +  2yhi  —  y2)0$/ 2  —  (hi  +  h2  —  y)0 7  +  h2(y  —  hx)09  -f  2x(hi  —  y)0io  + 
x0n  -xh20x 2  +  x20l3 

t xy  =  (y-  t*i)02  +  ( y 2  -  h\)0 4  +  (a  -  x)07  +  (x2  -  a2)0io  +  h2(a  -  x)09  +  x(y  -  Aj  )03 

a2Xx  -  a(x  ~  a)09  +  (x2  -  a2)0u  +  y(x  -  a)0Xi 

<*ly  =  0s  +  x0xi  +  y07  +  xy0x2  +  x20x3  (56) 

Txy  =  («  ~  x)07  -I-  ( x 2  -  a2)0xo  +  y(x  -  a)09 

ffxx  =  016  -  X0n  +  y0xs  -  xy0xg  -  x2i320/2  +  y202x 

aly  =  —(A3  +  2h3y  +  y2)02o/2  +  (h2  +  A3  +  y)0 7  +  0s  -  h2(h3  +  y)09  —  2x(A3  +  y)0 10  + 

x0ii  +  h2x0x2  +  x20x3 

Txy  =  (A3  +  y)0n  +  (y2  —  h\)0x9  +  (a  —  x)07  +  (x2  —  a2)0 10  +  h2(x  —  a)0 9  +  x(h3  +  y)02o 
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6.4  The  H2L13N  element 


The  H2L13N  element,  as  shown  in  Figure  50,  is  formulated  as  a  2-layered  element  incorporating  quadratic 
displacement  fields  resulting  in  8  nodes  per  layer.  Of  interest  in  this  element  formulation  is  the  inherent 
improvement  in  the  representation  of  in-plane  bending  behavior  due  to  the  quadratic  isoparametric  displace¬ 
ment  fields  in  each  layer  which  is  lacking  in  the  linear  order  elements  such  as  II2L6N  and  incompletely 
represented  in  H2L10N.  Incompatible  displacement  modes  can  be  added  to  the  linear  order  elements  to 
improve  bending  response  but  have  not  been  included  in  the  present  effort. 


Figure  50.  H2L13N  element  configuration. 

An  unexpected  feature  in  the  development  of  higher-order  layered  element  configurations  is  the  presence 
of  spurious  zero  energy  modes  in  the  resultant  element  stiffness  matrix  which  cannot  be  controlled  by  adopt¬ 
ing  higher-order  stress  expansions  while  enforcing  ail  layer  equilibrium  and  interface  continuity  conditions. 
In  order  to  qualitatively  assess  this  effect,  various  stress  fields  were  developed.  The  strain  field  obtained  from 
the  quadratic  isoparametric  displacement  field  yields  23  fundamental  strain  modes  which  interact  with  the 
assumed  stress  fields  in  generating  element  stiffness  characteristics.  In  following  the  procedure  for  selecting 
stress  expansions  in  the  linear-order  2-D  element,  complete  cubic  and  quartic  stress  fields  were  developed 
which  incorporate  28  and  40  independent  stress  modes,  respectively,  yet  both  assumed  fields  yield  spurious 
kinematic  modes  in  the  resulting  element  stiffness  matrix.  The  cubic  order  expansion  is  identical  to  that 
used  in  H2L6N  and  is  presented  in  equation  (43).  The  quartic  order  expansion  is  not  shown  due  to  the 
inadequacy  of  the  assumed  field  to  suppress  spurious  zero  energy  modes.  Although  the  individual  element 
stiffness  matrices  are  rank  deficient,  sufficient  constraints  are  introduced  through  the  assembly  and  sup¬ 
port  conditions  of  the  complete  model  to  contain  the  propagation  of  zero  energy  modes.  Predictions  using 
H2LI3N  incorporating  the  complete  cubic  order  stress  field  are  presented  in  Figures  51  and  52. 
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Figure  51.  II2L13N  tryy  stress  distributions  incorporating  complete  cubic  stress  expansions. 


Figure  52.  H2L13N  rry  stress  distributions  incorporating  complete  cubic  stress  expansions. 


As  shown,  the  solution  for  normal  peel  stress  distribution  is  accurate  but  the  transverse  shear  stress  pre¬ 
diction  is  shown  to  converge  slowly  to  the  reference  solution.  To  assess  the  effect  of  increasing  the  order  of 
assumed  stress  expansions,  predictions  using  H2L13N  incorporating  complete  quartic  order  stress  field  are 
presented  below  in  Figures  53  and  54. 
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In  an  attempt  to  improve  element  behavior,  various  augmentations  of  the  applied  stress  constraints  were 
assessed  and  are  detailed  below.  A  version  of  II2L13N  was  developed  in  which  continuity  constraints  were 
violated  in  selected  higher-order  stress  terms  to  obtain  the  requisite  number  of  rigid  body  modes.  A  detailed 
analytic  approach  as  presented  in  Reference  [32]  was  used  to  precisely  identify  the  strain  modes  which  are 
independent  of  the  constrained  stress  field  and  suggested  particular  higher-order  terms  to  be  added  to  the 
assumed  stress  expansions  in  order  to  suppress  zero  energy  modes.  Continuity  is  violated  by  the  addition  of 
two  terms  -  0i$x*  and  fin*4  -  to  the  shear  stress  expansions.  The  influence  of  higher-order  terms  have,  in 
general,  a  diminishing  influence  on  recovered  stresses  yet  function  to  stabilize  the  element  stiffness  matrix. 
The  resulting  stress  field  contains  27  independent  stress  modes  and  may  be  contrasted  with  the  complete 
cubic  field  presented  in  Equation  (43).  The  developed  stress  field  is  given  by 

=  01  +  02*  +  03y  +  04  X2  +  foxy  +  06y2  +  fax3  +  Z0gx2y  +  3  foxy2 
<Tyy  =  0n  +  0\(h\  —  2/»iy  +  y2)  +  0&(2h3  —  ih2y  +  y3)  +  0i6(y  —  hi  —  h2)  +  0nh2(h2  +  2/»i  —  2y)  + 

023 A2(3y  —  h2  —  3Ai)  +  Z07x{h\  —  2h\ y  +  y2)  +  02ox(y  —  hi  —  h 2)  -f  20iixhi(2h\  +  h2  —  2 y)  i- 

014*  +  02s3*2(y  —  hi  —  h2)  +  018*2 

Txy  =  012  +  02(/*i  —  y)  +  0s(h2  -  y2)/2  +  0g (h3  -  y3)  +  01^2  -  0igh2/2  +  /?24 h3  + 

2 04x(hi  -  y)  +  0gx(2h2  -  3 y2)  -  0i6x  +  20nh2x  -  Z023 h\x  +  307 x2(ht  -  y)  -  (57) 

02oX2/2  +  3022h2X2  —  /?25*3  +  026*4 

<T2j.  =  /3l0  +  013*  +  Asy  +  017*2  +  019*y  + /^2iy*  +  022*3  +  3/?23*‘y  +  3^24*!/“ 

<Xyy  =  011  +  014*  +  016y  +  0ny2  +  018*“  +  020*3  +  3/?22*y2  +  023 y3  +  3/?25*2y 

Try  =  012-013y-016*-2/?i7*y-019y2/2-/?2O*2/2-3^22*2y-3/?23*y2  - 
024y3  -  025*3  +  027*4 

Figures  55  and  56  depicts  normal  and  shear  stress  predictions  along  the  adhesive/adherend  interface  using 
H2L13N  incorporating  the  above  stress  field. 
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Figure  55.  H2L13N  <ryy  stress  distributions  with  continuity  violated  in  higher-order  shear  stress  terms. 
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Figure  56.  H2L13N  rxy  stress  distributions  with  continuity  violated  in  higher-order  shear  stress  terms. 


The  eigencharacteristics  of  the  element  demonstrate  the  requisite  three  rigid  body  modes  and  the  conver¬ 
gence  is  improved  for  both  stress  components.  The  peak  normal  and  shear  stress  predictions  at  the  joint 
end  -  which  is  a  nonconvergent  point  in  the  model  in  linear  elastic  analysis  -  shows  that  this  version  of 
H2L13N  is  superior  in  ‘capturing’  the  singular  behavior  than  even  the  highly  discretized  model  used  for  the 
reference  solution.  Although  strict  continuity  is  violated  in  the  interface  shear  stresses,  the  unconstrained 
shear  expansion  terms  are  of  quartic  order  which  have  only  a  small  contribution  to  the  total  recovered  shear 
stresses  yet  are  capable  of  stabilizing  the  element  by  removing  spurious  zero  energy  modes. 


In  order  to  avoid  violating  continuity  constraints,  an  alternative  formulation  was  developed  in  which  full 
interface  continuity  is  enforced  but  layer  equilibrium  conditions  are  relaxed  such  that  the  field  equilibrium 
equations  were  satisfied  in  an  integral  sense  rather  than  pointwise.  This  condition  is  enforced  by  requiring 
that 


+  —jdv  =  0, 
6y 


+  6^k]dv  =  o 


(58) 


With  this  formulation,  stress  field  constraints  due  to  enforcing  equilibrium  are  significantly  reduced.  The 
resulting  field  contains  26  independent  stress  modes  given  by 


viz  =  06  +  y20i  +  x20?  +  xy03  +  y0A  -  x0$ 

°VV  =  014  + (y2  —  hl)0T  +  {y  —  hl)0%  + *{y  —  hl)09  +  X201O  —  hlX0u  —  tl20i2  + 

2013  + 

Txy  —  019  +  x(y  ~  h\)/hi0s  +  y20\6  —  xy/hy0i2  +  x20\7  +  y0s  -  xyh^/hiPig 

viz  ~  020  +  y202l  +  X"  022  +  xy0 23  +  y024  +  x025  (59) 

v^y  =  0i4  +  x20io  +  xy0 1 1  +  y0i2  +  x0i3  +  y2/?is 

Txy  —  0\9  +  {y2  —  h\)02O  —  (h"i  + y)02S  + h\0\6  —  X0\2+ X20n  +  h\06  +  xy0\s 
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Figures  57  and  58  depicts  the  interface  stress  predictions  using  1I2L13N  incorporating  the  above  stress  field. 


0.0  1.0  2.0  3.0  4.0  5.0  6.0 
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Figure  57.  H2L13N  <ryy  stress  distributions  with  relaxed  equilibrium  constraints. 


Figure  58.  I12L13N  tx y  stress  distributions  with  relaxed  equilibrium  constraints. 
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An  eigenanalysis  of  this  dement  demonstrates  the  presence  of  two  spurious  zero  energy  modes  indicating 
that  the  relaxation  of  equilibrium  wi..  the  ‘weak’  form  as  given  above  in  equation  (59)  still  yields  stress 
fields  which  do  not  span  the  strain  space.  The  performance  of  the  element  is  good  in  the  particular  single-lap 
joint  used  to  provide  a  comparison  between  various  elements  assessed  in  the  current  study  but  is  unsuitable 
for  use  as  a  general  element. 

To  explore  in  greater  detail  the  effect  of  imposed  constraints,  a  version  of  H2L13N  was  developed  in  which  all 
interface  continuity  conditions  were  eliminated  and  only  layer  equilibrium  constraints  were  enforced  point- 
wise.  A  complete  quadratic  field  yielded  three  spurious  kinematic  modes  and  was  discarded.  A  complete 
cubic  field  was  found  to  span  the  strain  space  which  yielded  only  the  requisite  zero  energy  modes  associated 
with  rigid  body  modes.  The  resulting  field  contains  36  independent  stress  modes  given  by 

ff\*  -  0io  -  Zzy'01  +  y'0 2  +  ylh-  *2y04  -  *3Ps/$  -  2* yp6  -  x2p7/2  -  xpa  -  y3/?9 

=  017  -  3x2t/0n  -  y3 13^/3  -  xy2fa  -  2xyp\7  -  y2P7/2  -  yfi13  +  x3/?H  +  *20is  +  */?ie 
rj,  =  Pis  +  y30i  +x3  +  xy2P4  +  x2yPs  +  y2P6  +  x2p\2  +  xyp7  4-  y0s  +  xpl3 

=  p26-x2yPl9-x3p2o/3-2xyp7l+y2P2i-x2p-23/2-xP24  +  y02s-2xy2p27  +  y3P7S  (60) 

=  03S-  y0i9  -  Zx2yp30  -  y3Pis/3  -  xy2P20  +  x3p3l  -  2xyPz2  -  y2P 23/2  +  x2/?33  +  */?34 
r2v  =  P36  +  y3Pn  +  x3Ps  0  -f  xy2Pi9  +  x2yp2  0  +  y2P2 1  +  x2p32  +  xyP23  +  yp2A  +  xp^  9 

Predictions  of  bondline  stresses  are  depicted  below  in  Figures  59  and  60. 
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Figure  59.  H2L13N  <rsy  stress  distributions  with  only  layer  equilibrium  conditions  enforced. 
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Figure  60.  H2L13N  Txy  stress  distributions  with  only  layer  equilibrium  conditions  enforced. 


As  shown,  the  normal  stress  predictions  show  a  convergence  to  the  reference  solution,  but,  with  the  removal 
of  continuity  constraints,  the  shear  stress  distribution  is  incorrectly  predicted. 

An  element  utilizing  a  complete  quadratic  stress  field  was  formulated  in  which  the  enforcement  of  layer 
equilibrium  was  eliminated  and  only  interface  continuity  was  enforced.  As  presented  below,  the  resulting 
stress  field  possesses  30  independent  stress  modes  and  does  not  exhibit  any  spurious  zero  energy  deformation 
modes. 


=  06  +  y2 0\  +■  x202  +  xy0 3  +  y0A  4-  x05 
~  09  +  *(y  -  hi)07  +  y0a  +  x'0io  ~  xhiju  +  x0l2  +  y2/?i3 

Txy  =  019  +  y~  014  +  *y0l5  +  X~  016  +  y0\~  +  X0ls 

axx  —  020  +  y"  021  +  x~022  +  xy023  +  1/024  +  x025  (61) 

Vyy  =  00  +  (jr  —  +  ( /*  2  +  y)027  +  ll\0a  +  X-010  +  X1J0H  +  X012  +  h\013 

T\y  ~  019  +  (l/2  ~  hn)028  +  x(y  +  1*2)029  +  (/*2  +  *j)03O  +  h\0\4  +  h\x016  + 

x"  0i6  +  f*i0n  +  x0\s 


The  resulting  element  stress  predictions  are  presented  below  in  Figures  61  and  62. 
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Figure  61.  H2L13N  <ryy  stress  distributions  with  only  continuity  conditions  applied. 
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Figure  62.  H2L13N  rxy  stress  distributions  with  only  continuity  conditions  applied. 

Rapid  convergence  to  the  reference  solution  is  demonstrated  which  constrasts  with  the  reduced  conver¬ 
gence  behavior  demonstrated  in  H2L6N  with  similar  stress  field  constraints  applied.  This  shows  that  in  a 
higher-order  formulation,  the  neglect  of  satisfying  layer  equilibrium  appears  to  have  little  influence  in  accu¬ 
rately  predicting  stress  distributions  in  the  problem  under  study.  It  also  demonstrates  the  importance  of 
enforcing  stress  continuity  at  the  layer  interfaces  in  realizing  greater  accuracy  in  layered  hybrid  formulations. 
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A  comparative  eigenanalysis  is  presented  in  Table  2  below.  Isotropic  material  properties  were  selected  with 
E  =  1000.0  and  fi  =  0.25  for  each  layer  and  layer  dimensions  were  chosen  to  be  unity.  For  a  2-D  element,  the 
three  required  rigid  body  modes  are  depicted  in  Figure  63  and  the  two  spurious  kinematic  modes  resulting 
from  the  use  of  cubic  and  quartic  stress  fields  are  depicted  in  Figure  64.  Without  performing  a  detailed 
analysis  of  independent  element  strain  energy  contributions,  it  is  apparent  that  the  combined  application 
of  field  equilibrium  and  interface  continuity  conditions  constrain  the  assumed  stress  field  such  that  spurious 
kinematic  modes  are  unavoidable  without  selective  relaxation  or  violation  of  constraints  to  stabilize  the  ele¬ 
ment.  The  continuity  constraints  link  the  stress  fields  in  each  layer  and  tend  to  combine  stress  modes  which 
are  a  function  of  the  thickness  coordinate  which  diminishes  their  span  of  the  strain  space  as  a  function  of  y. 
As  can  be  seen  in  Figure  64,  the  deformation  modes  corresponding  to  zero  strain  energy  states  qualitatively 
show  that  they  are  exclusively  composed  of  v  component  deformations  which  are  functions  of  the  thickness 
coordinate. 


Table  2.  Comparison  of  eigenvalues  obtained  in  the  II2L13N  element  incorporating  various  stress 
fields:  complete  cubic  (Cubic-I),  complete  quartic  (Quart-I),  modified  cubic  (Cubic-II), 
quadratic  with  relaxed  equilibrium  constraints  (Quad-1),  cubic  with  only  equilibrium 
enforced  (Cubic  III)  and  quadratic  with  only  continuity  conditions  (Quad-II). 


1  A 

MwnffraB 

Quartic-I 

Cubic-II 

Quad- 1 

Cubic-Ill 

Quad-II 

1  1 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

23.69772 

0.0 

72.43028 

60.86598 

0.0 

0.0 

34.05380 

0.0 

163.0501 

78.74336 

6 

75.76981 

82.01182 

76.95349 

60.53644 

225.6264 

156.9013 

7 

126.6998 

227.4995 

134.2405 

81.16368 

268.1523 

161.8190 

8 

178.7381 

278.2503 

196.9742 

80.78365 

314.8713 

331.5214 

9 

211.6415 

337.8555 

218.3887 

172.7010 

388.8823 

358.9699 

10 

293.2313 

378.8533 

316.7293 

173.2961 

391.7162 

4 13. 1536 

11 

325.1373 

441.0916 

341.2829 

331.0209 

458.1095 

557.5755 

12 

355.5333 

480.2937 

386.6343 

376.4624 

553.0599 

564.5969 

13 

371.5992 

545.4102 

431.6523 

424.4714 

603.3747 

809.1861 

14 

405.7435 

568.8378 

530.7245 

474.6536 

733.2789 

920.7719 

15 

517.5139 

602.9718 

587.1902 

559.7191 

753.5590 

988.4729 

16 

587.8730 

715.4142 

640.0043 

577.7691 

766.3474 

1010.934 

17 

674.6511 

773.2681 

768.9178 

929.2014 

1040.206 

1259.228 

18 

766.4438 

152.7025 

901.8028 

913.4121 

1152.900 

1403.390 

19 

819.4111 

950.2543 

1013.777 

1071.387 

1221.521 

1643.767 

20 

1044.450 

1110.563 

1083.736 

1211.052 

1361.066 

1974.218 

21 

1087.227 

1181.065 

1123.716 

1389.199 

1420.998 

2549.088 

22 

1210.914 

1288.766 

1555.517 

1672.933 

1537.543 

2768.760 

23 

1469.354 

1548.891 

1575.495 

1726.760 

1851.332 

3645.037 

24 

1529.672 

1595.397 

1580.724 

2022.144 

1969.350 

4235.017 

25 

1866.629 

1871.165 

1936.631 

2680.260 

2462.195 

4727.452 

26 

1950.449 

1961.606 

2113.241 

3173.282 

3038.038 

7906.456 
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TRANSLATIONAL  MODE  I  TRANSLATIONAL  MODE  II 


Figure  63.  Required  translational  and  rotational  rigid  body  modes 


ZERO  ENERGY  MODE  I  ZERO  ENERGY  MODE  II 

Figure  64.  Spurious  kinematic  modes. 


The  inability  to  enforce  all  stress  held  constraints  without  giving  rise  to  zero  energy  modes  makes  the 
general  utility  of  the  higher-order  H2L13N  element  questionable.  Although  accurate  stress  predictions  are 
demonstrated  in  the  ‘Cubic-II’  version  in  which  higher-order  expansion  terms  are  introduced  which  vio¬ 
late  continuity  and  in  the  ‘Quad-II’  version  in  which  only  continuity  is  enforced,  element  performance  may 
degrade  in  the  modelling  of  other  general  adhesive  joint  configurations. 

6.5  Remarks  on  2-D  Element  Behavior 

The  detailed  study  undertaken  with  various  2-D  element  formulations  has  highlighted  certain  features  of  the 
hybrid  stress  method  in  layered  element  configurations  which  may  provide  insight  into  the  use  of  the  hybrid 
technique  in  other  similar  applications.  For  the  present  use  in  developing  layered  element  formulations 
for  the  analysis  of  adhesive  joints,  various  observations  may  be  summarized.  The  study  of  the  linear- 
order  H2L6N  element  involving  selective  enforcement  of  stress  equilibrium  and  continuity  constraints  has 
demonstrated  that  the  best  performing  element  formulation  is  one  which  strictly  enforces  all  constraints. 
A  related  conclusion  is  drawn  from  the  use  of  the  partial  hybrid  stress  functional  as  the  stress  prediction 
has  been  shown  to  deteriorate  due  to  the  mixed  nature  of  the  element  incorporating  both  displacement- 
based  and  stress-based  element  field  assumptions.  Related  to  the  H2L6N  element  performance  in  which  held 
equilibrium  conditions  are  neglected,  the  inability  of  enforcing  held  equilibrium  in  each  layer  with  the  partial 
hybrid  functional,  leads  to  inaccuracies  in  the  prediction  of  bondline  stresses.  The  investigation  of  selective 
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higher-order  element  formulations  such  as  in  H2L10N  results  in  an  element  which  offers  no  advantage.  The 
formulation  of  a  fully  higher-order  layered  element  configuration  as  in  H2L13N  results  in  an  element  with 
an  inherent  presence  of  kinematic  modes  if  all  constraint  equations  are  strictly  applied  to  the  assumed 
stress  fields.  The  performance  of  various  higher-order  element  formulations  over  lower-order  elements  have 
not  demonstrated  significantly  improved  performance  in  the  prediction  of  bondline  stresses.  It  is  therefore 
concluded  that  the  optimum  2-D  stress  prediction  is  obtained  in  linear-order  layered  element  formulations 
in  which  all  stress  constraints  are  identically  satisfied  pointwise. 

7  3-D  Special  Adhesive  Elements 

Two  3-D  adhesive  joint  problems  are  analyzed  involving  a  rectangular  and  a  tapered  single-lap  joint.  Nodal 
variables  are  restricted  to  a,  v  and  w  translational  degrees  of  freedom.  The  joint  geometries  are  depicted 
above  in  Figures  5  and  6.  The  applied  loading  and  boundary  conditions  are  shown  above  in  Figure  7. 
As  previously  discussed,  reference  solutions  are  obtained  from  a  highly  discretized  2-D  model  using  8-node 
plane-strain  elements  shown  in  Figure  4.  For  the  3-D  single-lap  joint  with  rectangular  planform,  the  2-D 
reference  solution  is  used  directly  assuming  negligible  three-dimensional  effects  through  the  joint  width.  In 
the  3-D  tapered  joint,  stresses  from  the  2-D  reference  solution  are  linearly  scaled  to  account  for  the  effect  of 
stress  intensification  as  a  function  of  joint  taper.  The  reference  solutions  for  the  <rC2  and  tz1  stress  distribu¬ 
tions  together  with  a  depiction  of  the  Ty2  stress  field  over  the  adhesive  domain  are  presented  in  Figures  65 
through  70  for  the  rectangular  and  tapered  joint  configurations.  The  ryz  distribution  was  obtained  using  the 
developed  3-D  hybrid  elements  discussed  below  and  shows  the  three-dimensional  singular  nature  of  this  stress 
component  at  the  bond  corners.  Because  all  stress  components  are  related  through  equilibrium  constraints, 
the  presence  of  this  effect  causes  a  slight  departure  of  predicted  arz2  and  rxi  stresses  from  the  reference 
solution  obtained  from  the  2-D  plane-strain  model.  Stress  recovery  for  the  rectangular  configuration  was 
obtained  along  the  lower  adhesive/adherend  interface  while,  in  the  tapered  joint,  stresses  were  recovered 
along  the  upper  interface. 
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Figure  65.  Reference  atl  stress  distribution  in  the  rectangular  3-D  single-lap  joint. 
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Figure  66.  Reference  rgx  stress  distribution  in  the  rectangular  3-D  single-lap  joint. 


SINGLE-LAP  JOINT  (RECT.) 
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Figure  69.  Reference  rxt  stress  distribution  in  the  tapered  3-D  single-lap  joint. 


Figure  70.  Singular  r„,  stress  field  in  the  tapered  3-D  single-lap  joint. 

From  the  detailed  study  of  various  2-D  element  performance,  the  optimum  element  configuration  for  use  in 
general  3-D  applications  is  selected  as  a  two-layer  formulation  incorporating  bilinear  isoparametric  assumed 
displacement  fields  and  using  the  complete  Ilellinger-Reissner  functional  in  which  all  field  equilibrium  and 
interface  continuity  conditions  are  identically  satisfied.  No  convergence  studies  are  performed  in  the  3-D 
case;  the  degree  of  discretization  is  fixed  at  100  elements  along  the  bondline  with  5  elements  through  the 
width.  The  lack  of  mesh  refinement  through  the  joint  width  is  due  to  the  assumption  of  negligible  3-D 
stress  variations.  Results  are  compared  with  stress  predictions  using  8-node  C3D8  displacement-based  solid 
elements  available  from  the  ABAQUS  library  and  the  reference  solution. 


7.1  The  H2L12N  Element 

The  H2L12N  element  is  formulated  as  a  12-node,  two  layer  element  as  shown  in  Figure  71. 


Figure  71.  The  H2L12N  element  configuration. 
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As  a  3-D  extension  of  the  H2L6N  element,  two  version  j  were  developed  in  which  complete  linear  and  quadratic 
order  expansions  were  assumed  for  each  stress  component.  Equilibrium  is  enforced  within  each  layer  and 
continuity  of  the  a„,  txi  and  ryj  stress  components  imposed  at  the  layer  interface.  The  H2L12N  element  is 
formulated  to  represent  a  general  planform  geometry  to  provide  a  robust  element  for  3-D  analysis  of  arbitrary 
adhesive  joint  configurations. 

The  linear  stress  field  contains  48  independent  stress  modes  and  is  given  by 

*1,  =01-  *03  -*03  +  y04  +  Z0S-  *y09  -  **07  +  *V0a 
=  09  +  *010  +  y0u  +  * 0i3  +  *y0i3  +  **0i4  —  *y0n 
9\i  —  0ia  +  (hi  +  hi  —  z)0u  +  (hi+h2  —  g)0ir  +  t(hi  +  h3~’  z)0i9  +  z02o  + 
y(hi  +  hj  —  z)0j  i  +  y023  +  *y03a 

Tyz  =  037 +  (*  — hl)02i  —  1*3038  + *(hl  -  *)013  + *029  +  fl2*03O  +  y0n  +  *y0l9 

*zx  =  031  +  (*  ~  hi)03  -  ho033  +  *016  +  y(*  ~  f*l)06  +  V033  ~  f*2y034  +  *y02l  (62) 

Tly  =  024  ~  * 011  ~  * 02 S  +  V03  +  *026  +  **016  +  *y07 

9lt  =  036-*036  —  *032  +  y037  +  *038-*y034-**039  +  y*04O 

<7yy  =  041  +  *042  +  y047  +  *043  +  *V03O  +  **044  —  *V046 

9z,  =  018  + *020  + V022  ~  *016  ~  *017  +  *y023  ~  **019  ~  *y021 

Tyt  =  027  +  *029  +  V017  +  *028  +  *V019  ~  **030 

~  031  +  *016  +  y033  +  *0 32  +  Xy021  +  y*034 

*xy  =  046  ~  * 047  ~  *028  +  V036  +  * 048  +  **046  +  *y039 

The  quadratic  stress  field  contains  78  independent  stress  modes  and  is  given  by 

9ix  =  010  +  y*01  +  *203  +203+  y204  +  X~06  +  *V06  +  V07  +  *08  +  *2  09 
9W  =  0l9  +  y*01i+*Z012  +  Z013  +  y20l4  + X1013  +  *y016  +  y0n  + X018  + *202O 

=  030  +  *(*  —  2Al  )021  +  yz0 22  +  XZ0 23  +  Z024  +  y2 02 5  +  X2 026  +  *y027  +  V0 28  + 

*029  —  2zh3031 

Ty,  =  043  +  »(2Al  -  Z)0 21  ~  *V023  ~  V024  ~  y*014  ~  1x*034  ~  **016  ~  *0 17  -  V041  + 

*042  -  2036  +  2yh2031  +  Z20 39  +  y*06  +  y2 044  -  2xy043  +  X2 046 
Tlz  =  048  ~  *y022  ~  %yZ033  +  XZ0 M  +  y0A7  -  Z0 35  +  X0A1  +  Z2  038  ~  **{■ ?5  -  yZ06  - 
XZ021  ~  Z08  +  y2049  ~  "l*y044  +  X2 046 

T^y  =  037-y*02-XZ011+Z032  +  y2033-Xy014+X2034+y03S  +  X036-2yZ038-  (63) 

2xz039  -  xy0i  +  z2 04o  +  xy02i 

9xz  =  068  +  y*06O  +  XZ0 51  +  Z032  —  2xy0$3  +  y2  064  +  X2066  +  y0 56  +  1067  +  Z2 0 59  -  Xy0 60 
9yy  =  069  +  yz0 61  +  XZ0 62  +  *063  +  y2 064  +  X2 086  +  *y066  +  V067  +  *068  +  *2 070 

<f2t  =  030  +  (~h|  +  z2  —  2hih2)03i  +  y(/»2  +  hi  +  z)02 2  +  x(h 2  +  hi  +  z)0 23  -I-  (hj  +  hi  +  z)0 24  + 

y2028  +  X2026  +  xy027  +  V028  +  *029  ~  h202 1 

Ty,  =  043  +  (z2  -  h\)0 77  +  y(z  +  h2)0 55  —  (^2  +  *)067  +  y(*2  ~  *)031  ~  y{*  +  h2)064  ~  x(z  +  *2)^66  - 

(/»2  +  z)07 3  —  2x(z  +  h2)072  —  *V023  ~  V0 24  ~  hiy0\4  —  2hxX034  —  hiX0 16  —  hi0n  —  y04i  + 

*042  —  hl038  +  h2039  +  h\y06  +  hiy021  +  y2044  —  1*y048  +  X2046 
*iz  —  048  +  (z2  -  h\)028  +  y(z  +  *2)^60  +  *(*  +  h2)064  —  *(z  +  ^2)^55  -  l(z  +  *2^31  +  (h2  +  z)0 75  - 

*y022  -  2hiy033  +  hix0!4  +  y0A7  -  hi033  +  x0Ai  +  h.]03s  -  hix0h  -  hiy0e  -  hix02i  - 
hl08  +  y2049  —  2*y044  +  *2043 

T2y  =  074  -  *2061  -  y* 061  +  *071  +  y20 53  ~  Xy064  +  *2072  ~  *V066  +  X073  ~  V067  +  *V031  - 
y073  -  2y* 076  ~  2XZ077  +  Z2 078 
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The  effect  of  high  aspect  ratio  is  examined  using  the  element  configuration  shown  in  Figure  72.  Isotropic 
material  properties  were  selected  with  E  =  1000.0  and  /i  =  0.25.  An  eigenanalysis  of  the  two  versions  of 
H2L12N  compared  with  a  displacement-based  formulation  is  depicted  in  Table  3  and  shows  the  elimination 
of  low-energy  or  ‘weak’  characteristic  deformation  modes  by  incorporating  a  higher-order  expansion  for  the 
stresses.  In  addition,  the  hybrid  formulation  has  the  desirable  feature  of  relieving  the  stiff  deformation  modes 
demonstrated  in  the  displacement-based  element. 


Figure  72.  The  H2L12N  element  with  a  high  aspect  ratio  layer. 


Table  3.  Comparison  of  eigenvalues  obtained  in  the  II2L12N  element  incorporating 
complete  linear  and  quadratic  stress  expansions  with  ”  displacement-based 
element. 


Figures  73  through  78  depict  stress  predictions  over  the  bond  surface  and  along  the  centerline  for  a  3-D 
single-lap  adhesive  joint  with  a  rectangular  planform.  Solutions  are  obtained  for  the  normal  peel  stress 
distribution,  <7„,  and  the  transverse  shear  stress  distribution,  rzz  ,  across  the  bond  surface.  Quadratic  stress 
fields  are  assumed  in  the  1I2L12N  element  and  stress  predictions  are  compared  to  results  obtained  using 
displacement-based  elements  and  the  reference  solution. 
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Figure  73.  H2L12N  <rtz  stress  prediction  over  the  bond  surface  in  rectangular  3-D  single-lap  joint. 


Figure  74.  C3D8  (TZ1  stress  prediction  over  the  bond  surface  in  rectangular  3-D  single-lap  joint 
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SINGLE-LAP  JOINT  (RECT.) 


Figure  77.  C3D8  rxt  stress  prediction  over  the  bond  surface  in  rectangular  3-D  single-lap  joint. 


Figure  78.  rxz  stresses  along  centerline  of  adhesive  bond  in  rectangular  3-D  single-lap  joint. 

As  shown,  the  H2L12N  element  predictions  are  excellent  while  the  displacement-based  element  solutions 
converge  away  from  the  reference  solution.  The  slight  departure  from  the  reference  solution  is  attributed  to 
the  deviation  from  a  state  of  plane-strain  due  to  the  three-dimensional  nature  of  the  singularities  in  the  rvz 
stress  field  which  influences  crlz  and  rxz  through  the  equilibrium  constraints. 

Stress  predictions  are  obtained  for  the  normal  peel  stress  distribution,  <tzz,  and  the  transverse  shear  stress 
distribution,  rxz,  across  the  bond  surface  in  a  tapered  joint  and  are  shown  below  in  Figures  79  through  84 
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Figure  80.  C3D8  <rlt  stress  prediction  over  the  bond  surface  in  a  tapered  3-D  single-lap  joint. 
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Figure  82.  H2L12N  rrt  stress  prediction  over  the  bond  surface  in  a  tapered  3-D  single-iap  joint. 
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e  83.  C3D8  rX2  stress  prediction  over  the  bond  surface  in  a  tapered  3-D  single-lap  joint. 


Figure  84.  tX2  stresses  along  centerline  of  adhesive  bond  in  a  tapered  3-D  single-lap  joint. 
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7  2  Remarks  on  3-D  Element  Performance 

The  H2L12N  element  demonstrates  the  same  improvement  over  displacement-based  8-node  solid  elements 
as  was  shown  in  the  assessment  of  the  2-D  H2L6N  element.  The  3-D  element  formulation  developed  herein 
permits  a  general  adhesive  joint  planform  to  be  modelled  and  demonstrates  ideal  flexible  behavior  in  high 
aspect  ratio  element  configurations. 


8  Closure 

Numerous  specialized  layered  element  formulations  have  been  derived  to  fully  assess  the  performance  of 
the  hybrid  element  technique  in  the  application  to  improve  the  computational  efficiency  and  accuracy  of 
determining  critical  stresses  along  the  bond  interfaces  in  adhesive  joints.  The  best  performing  elements  have 
been  the  lower-order  element  formulations  due  to  the  difficulty  encountered  with  the  zero  energy  modes 
present  in  the  higher-order  layered  elements.  Optimal  element  performance  has  been  associated  with  the 
strict  enforcement  of  all  layer  equilibrium  and  interface  continuity  constraints.  A  clear  improvement  in  stress 
prediction  has  been  demonstrated  over  the  use  of  displacement-based  element  which  validates  the  superior 
behavior  of  the  hybrid  stress  technique  in  this  particular  application.  For  a  robust  finite  element-based 
analytical  framework  for  the  study  of  adhesive  bond  stresses,  the  optimum  element  configurations  developed 
herein  are  expected  to  form  the  basis  for  further  development.  Future  work  is  suggested  in  the  incorpora¬ 
tion  of  incompatible  displacement  modes  in  2-D  lower-order  elements  to  enhance  the  derived  elements  for 
problems  involving  significant  bending  stresses.  Nonlinear  material  behavior  should  be  supported  in  future 
efforts  to  predict  the  reduction  of  peak  stress  due  to  plastic  material  yielding  at  the  joint  ends  together 
with  the  capability  to  mode!  large  displacements  which  occur  in  various  adhesive  joint  configurations  under 
peak  service  loading.  Finally,  a  local  failure  analysis  methodology  should  be  developed,  perhaps  including  a 
micromechanical  level  analysis  of  the  adhesive/adherend  interface,  to  predict  the  initiation  of  bond  failure. 
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